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PREFACE 


This  document  was  written  to  serve  as  an  addendum  to  Master's 
Thesis  project  AFIT/GNE/ENP/91M-6,  "High  Altitude  Neutral 
Particle  Transport  Using  the  Monte  Carlo  Simulation  Code  MCNP 
With  Variable  Density  Atmosphere".  The  research  and  develop¬ 
ment  for  this  thesis  was  provided  by  Captain  David  L.  Monti, 
USAF.  Assistance  was  provided  by  Lieutenant  Commander  Kirk 
A.  Mathews,  USN.  This  document  was  written  to  provide  under¬ 
standing  and  assistance  to  users  who  plan  to  implement 
additional  modifications  to  MCNP. 

This  effort  was  part  of  an  ongoing  AFIT  research  project  to 
imp-'-ovu  accuracy,  decrease  run  time,  and  improve  overall 
Monte  Carlo  transport  methods  at  AFIT.  There  were  two  primary 
goals  of  this  thesis  project.  The  first  was  to  develop  an 
accurate  model  describing  the  variation  of  air  density  with 
aitituoe  and  to  incorporate  it  into  a  generalized  version  of 
tne  Monte  Carlo  code  MCNP.  The  second  goal  was  to  perform 
r'aiiation  transport  simulations  using  a  point  isotropic 
neutron-photon  source  and  study  the  effects  over  long  ranges. 

This  work  was  sponsored  by  the  Air  Force  Weapons  Laboratory; 
this  AFIT  technical  report  is  intended  to  serve  as  a  formal 
addendum  to  the  actual  thesis  project.  A  primary  goal  in 
this  report  is  to  provide  the  sponsor  and  all  future  users  of 
the  modified  MCNP  code  with  detailed  information  regarding 
the  additions  and  modifications  made  to  MCNP,  version  3E.  We 
hope  it  provides  the  necessary  documentation  to  assist  in 
future  work  in  this  area. 
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ABSTRACT 


^rate  a  continuously 


MCNP  version  3B  was  modified  to  incorp^ 
variable  density  atmosphere.  This  was,  accomplished  by  repre¬ 
senting  the  variation  of  air  density  ajs  a  -function  of 
altitude  with  a  series  of  continuous  pfiecewise  exponential 
curves  up  to  a  maximum  altitude  of  100|o  km.  User — written 
subroutines  and  functions  were  written  which  incorporated 
these  piecewise  functions.  These  subroutines  and  functions 
were  subsequently  incorporated  into  a/production  version  of 
MCNP.  Several  MCNP  subroutines  and  files  were  modified  in 
support  of  these  modifications.  This/  report  discusses 
detailed  information  regarding  the  theoretical  development  of 
the  variable  density  model,  the  userf-wri tten  subroutines  and 
functions,  the  modifications  to  MCNP  subroutines  and  files, 
and  other  relevant  information. 
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1  Introduction 


MCNP  version  3B  Cl  3  was  modified  to  incorporate  a  continu¬ 
ously  variable  density  atmosphere  as  part  of  a  Master’s 
Thesis  project.  This  addendum  is  part  of  Master's  Thesis 
AFI T /GNE/ENP/9 1 M-6  C23,  still  unpublished  at  this  writing. 
There  were  several  reasons  for  incorporating  these  modifica¬ 
tions,  These  include!  1)  improving  the  accuracy  of  computed 
Monte  Carlo  density-dependent  radiation  transport  results 
within  the  atmosphere,  2)  decreasing  the  overall  computer  run 
time,  3)  increasing  the  dimensions  of  the  problem  geometry, 
and  4)  lessening  the  user  overhead  in  preparing  the  input 
files.  A  full  discussion  on  this  topic  can  be  found  in 
Master's  Thesis  AFIT/GNE/ENP/91M-6  by  David  L.  Monti. 

There  are  several  areas  requiring  further  work  on  the  modi¬ 
fied  MCNP  code,  some  of  which  were  discussed  in  Section  VIII 
of  AFIT/GNE/ENP/91M-6.  This  addendum  was  written  to  assist 
future  users  in  performing  additions  and/or  modifications  to 
the  modified  MCNP  version  3B  code  at  the  Air  Force  Institute 
of  Technology  (AFIT). 

1.1  Background 

The  purpose  of  the  thesis  project  was  to  perform  Monte  Carlo 
simulations  of  neutral  particle  transport  with  primary  and 
secondary  photon  production  as  it  applies  to  a  point  iso¬ 
tropic  source  within  a  variable  density  atmosphere.  Specifi¬ 
cally,  the  goal  was  to  increase  the  accuracy  of  neutral 
particle  transport  calculations  within  the  atmosphere  by 
accurately  representing  the  variability  of  air  density  with 
altitude.  In  previous  studies  of  Monte  Carlo  radiation 
transport  within  the  atmosphere,  many  discrete  spatial  cells 
were  used  to  represent  the  change  in  air  density  as  a  func¬ 
tion  of  altitude.  This  was  accomplished  by  using  a  series  of 
constant  density  regions,  usually  with  a  constant  scale 
height  factor.  This  prol i f eration  of  cells  proved  to  be 
costly  in  terms  of  computer  run  time,  especially  for  large 
dimension  problems. 

During  this  study,  analytic  representations  for  air  density 
and  mass  integral  which  modelled  the  variation  of  air  density 
with  altitude  were  developed.  The  generalized  Monte  Carlo 
code  MCNP,  version  3B,  was  modified  using  the  newly  developed 
atmospheric  model.  The  work  included  the  development  of 
user-written  subroutines  and  functions  as  well  as  modifi¬ 
cation  of  MCNP  subroutines  and  files. 

1.2  Approach  to  Problem 

The  approach  to  the  problem  of  developing  a  variable  density 
atmosphere  was  based  on  a  two-phase  approach:  1)  develop 
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analytic  expressions  to  represent  both  the  variation  of  air 
density  with  altitude  and  cumulative  mass-integral  with 
altitude,  and  2)  implement  these  analytic  functions  into  MCNP 
via  user — written  subroutines  and  functions.  It  was  also 
necessary  to  perform  auxiliary  modifications  to  MCNP  files 
and  subroutines  in  support  of  the  variable  density  model. 

2  Project  Goals 

The  overall  objectives  of  the  thesis  project  were  as  follows: 

1.  Develop  a  working  Monte  Carlo  simulation  code  (MCNP)  that 
incorporates  a  continuously  variable  air  density  capability; 

2.  Improve  the  accuracy  of  density-dependent  transport 
results  over  previous  results; 

3.  Decrease  computer  run  time  relative  to  the  unmodified 
Monte  Carlo  code,  MCNP; 

4.  Increase  the  capability  to  run  very  large  dimension 
problems; 

5.  Study  the  radiation  effects  from  a  high  altitude  burst 
over  long  ranges. 

The  overall  objectives  of  this  report  are  as  follows: 

1.  Provide  detailed  information  on  the  atmospheric  density 
model ; 

2.  Provide  detailed  information  on  the  user-written  subrou¬ 
tines  and  functions  incorporated  into  MCNP; 

Provide  detailed  information  on  MCNP  files  and  subroutines 
modified  in  support  of  the  variable  density  atmosphere. 

3  Sequence  of  Presentation 

Section  4  presents  the  relevant  MCNP  program  flow.  Section  5 
provides  the  theoretical  discussion  behind  the  variable 
density  representation  of  the  atmosphere.  Section  6 
presents  a  detailed  discussion  of  the  user-written  subrou¬ 
tines  and  functions  integrated  into  MCNP  as  part  of  the 
modifications.  Section  7  provides  a  detailed  discussion  of 
the  affected  MCNP  subroutines.  Section  8  presents  some 
miscellaneous  information  concerning  the  MCNP  code.  Appendix 
A  contains  listings  of  the  user — written  subroutines  and 
functions  integrated  into  MCNP.  Appendix  B  contains  listings 
of  the  MCNP  files  and  subroutines  that  were  modified. 
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4  MCNP  Program  Flow 


This  section  describes  the  relevant  portions  of  MCNP  program 
flow.  Topics  discussed  include  the  program  flow  of  Subrou¬ 
tines  HSTORY  and  TRANSM,  where  most  of  the  particle  transport 
takes  place. 

4.1  Transport  Program  Flow 

Particle  transport  in  MCNP  is  controlled  via  the  main  trans¬ 
port  overlay  MCRUN.  MCRUN  calls  subroutine  TRNSPT  to  start 
particle  histories  via  subroutine  HSTORY. 

4.1.1  Subroutine  HSTORY 

Subroutine  HSTORY  begins  the  particle  history  by  starting  a 
particle  from  the  source  via  subroutine  STARTP.  All  source 
parameters  such  as  initial  direction,  weight,  and  energy  are 
determined.  Initial  source  parameters  are  next  checked 
against  user — supplied  weight  cut-offs  and  energy  cut-offs  and 
some  summary  information  is  incremented.  Subroutine  TALLY  is 
called  to  score  tally  contributions  during  a  particle 
history.  If  point  detectors  or  ring  detectors  are  being 
used,  then  subroutine  TALLYD  is  called  from  TALLY  to  score 
contributions  to  detectors.  If  DXTRAN  spheres  are  being 
used,  then  DXTRAN  is  called  to  score  DXTRAN  sphere  contribu¬ 
tions.  The  weight  window  variance  reduction  game  is  played 
if  requested,  and  energy  splitting  or  Russian  roulette  is 
performed  if  required. 

Back  in  subroutine  HSTORY,  particle  transport  continues  by 
calling  TRACK  to  determine  the  intersection  of  the  particle 
trajectory  with  each  bounding  surface  of  the  cell.  The 
distances  to  the  nearest  cell  boundary,  DLS,  time  cut-off, 
DTC,  and  DXTRAN  sphere,  DXL,  are  calculated.  Cross  sections 
for  the  cell  are  determined  via  subroutines  ACETOT  for 
neutrons  and  PHQTOT  for  photons.  A  macroscopic  cross 
section,  QPL,  is  calculated  and  then  modified  by  the  expo¬ 
nential  transform  in  subroutine  EXTRAN,  if  necessary.  The 
distance  to  next  collision,  PMP,  is  determined  by  a  random 
number  sampling  of  a  logarithmic  distribution.  Up  to  this 
point,  no  modif ications  were  made  to  subroutine  HSTORY. 

The  track  length  in  the  cell  is  determined  by  taking  the 
minimum  of  the  distance  to  collision,  PMF,  the  distance  to 
time  cut-off,  DTC,  the  distance  to  the  nearest  bounding 
surface,  DLS,  and  the  distance  to  a  DXTRAN  sphere,  DXL.  The 
particle  will  undergo  a  collision  only  if  the  distance  to 
collision,  PMF,  is  less  than  the  distance  to  a  surface  cros¬ 
sing,  DLS.  Subroutine  TALLY  is  called  to  increment  any  cell 
tallies  if  necessary  and  some  summary  information  is  incre¬ 
mented.  The  particle's  position  is  also  updated.  The 
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energies  and  directions  of  the  particles  exiting  a  collision 
site  are  handled  in  subroutines  ACECAS  and  ACEOS. 

4.1.2  Subroutine  TRANSM 

Subroutine  TRANSM  is  called  •from  TALLYD  to  calculate  the 
total  transmission  to  the  detector.  Prior  to  calculating  the 
transmission,  subroutine  DDDET  is  called  to  calculate  the 
distance,  DD,  to  the  detector  and  TRACK  is  called  to  calc¬ 
ulate  the  distance  to  the  nearest  intersecting  surface,  DLS. 

Cro=s  sections  for  the  cell  are  determined  via  subroutines 
ACETOT  for  neutrons  and  PHOTOT  for  photons.  Once  the  micro¬ 
scopic  cross  section  in  the  cell  has  been  determined,  it  is 
multiplied  by  the  atom  density  in  the  cell  to  obtain  the 
macroscopic  cross  section.  The  track  length  in  the  cell,  D, 
is  determined  by  taking  the  minimum  of  the  distance  to  the 
nearest  bounding  surface,  DLS,  and  the  distance  to  the 
detector,  DD.  The  number  of  mean  free  paths  to  the  detector 
or  nearest  bounding  surface,  AMFP,  is  determined  using  the 
values  of  the  macroscopic  cross  section  in  the  cell,  PLE,  and 
the  minimum  distance  to  the  detector  or  nearest  bounding 
surface,  D.  Then  the  locations  of  the  new  cell  and  particle 
position  are  updated. 

5  Theoretical  Discussion  of  the  Variable  Density  Atmosphere 

This  section  presents  a  detailed  theoretical  discussion  of 
the  variable  density  atmosphere  incorporated  into  MCNPi 
version  3B.  This  served  as  a  basis  for  the  development  of 
user-written  subroutines  and  functions  which  were  subse¬ 
quently  incorporated  into  MCNP.  Topics  discussed  include: 
definition  of  a  reference  density  and  calculation  of  impor — 
tant  atmospheric  transport  parameters. 

5.1  Reference  Density 

Modelling  the  variable  density  atmosphere  involved  defining  a 
reference  air  density.  In  the  unmodified  version,  MCNP 
computes  density-dependent  parameters  using  a  discrete  cell 
density  defined  in  the  input  file.  In  the  modified  version, 
the  variation  of  air  density  is  continuous,  and  a  constant 
density  spatial  mesh  in  the  vertical  direction  is  no  longer 
needed,  except  for  splitting/Russian  roulette.  Therefore,  in 
the  modified  version  of  MCNP,  the  air  density  in  each  cell  is 
redefined  to  be  the  mass  density  at  sea- level,  l.225xl0~3 
q/cm3,  regardless  of  the  values  read  from  the  input  file. 

All  subsequent  calculations  of  density-dependent  parameters 
will  now  be  calculated  based  an  a  sea-level  homogeneous 
atmosphere.  At  this  point,  prior  to  computing  the  distance 
to  collision  or  the  transmission  to  the  detector  or  nearest 
intersecting  boundary,  the  MCNP  calculations  are  intercepted. 
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The  aforementioned  parameters  are  then  corrected  for  a 
variable  density  atmosphere. 

Prior  to  running  the  MCRUN  overlay,  the  main  overlay  IMCN  is 
run  to  read  the  input  file  to  obtain  the  necessary  parameters 
used  in  defining  the  problem.  Subroutine  NEXTIT  is  called  to 
process  items  from  the  input  file,  i.e.,  to  store  such  para¬ 
meters  as  material  density  and  cell  descriptions.  Subroutine 
NEXTIT  was  modified  for  the  purpose  of  redefining  the  density 
in  each  cell  to  be  the  density  at  sea-level. 

5.2  Transport  Parameters 

In  Monte  Carlo  transport  calculations,  there  are  two  prin¬ 
ciple  parameters  that  involve  determining  the  optical  path 
length  between  two  points  in  space:  1 )  the  sampled  distance 
to  collision  and  2)  he  number  of  mean  free  paths  to  a 
detector  or  nearest  intersecting  boundary.  Both  parameters 
are  closely  related  and  depend  strongly  on  the  atmospheric 
density.  Details  of  the  calculation  of  each  of  these  two 
important  transport  parameters  is  presented  below. 

5.2.1  Distance  To  Collision 

The  distance  to  collision,  PMF,  computed  in  HSTORY,  is  based 
on  the  microscopic  cross  section  in  the  cell  and  the  atom 
density  in  the  cell,  which  was  redefined  in  NEXTIT  to  be  the 
air  density  at  sea-level.  This  is  given  by 

PLE-  TOTH*  DEN  ( CELL)  (1) 


where 


PLE  =  macroscopic  cross  section  in  current  cell  Ccm~‘3 
TQTM  =  microscopic  cross  section  in  current  cell  for 

either  a  neutron  or  photon  reaction  Cbarns/atom) 
DEN(CELL)  =  atom  density  in  current  cell  C a toms /barn -cm3 

Since  the  distance  to  collision  is  a  probability  based  on 
density,  it  is  necessary  to  correct  this  value  for  a  variable 
density  atmosphere.  Subroutine  EQDIST,  called  from  HSTORY, 
was  written  to  calculate  the  corrected  distance  to  collision 
based  on  values  input  from  HSTORY,  i.e.,  current  p.  tide 
altitude,  2-direction  cosine,  and  distance  to  collision  based 
on  a  sea-level  homogeneous  atmosphere.  Calculation  of  the 
required  quantities  used  to  determine  the  corrected  distance 
to  collision  are  explained  below.  See  Figure  3,  Section  III, 
AFIT /SNE/ENP/91M-6 ,  “Mass-Integral  Scaling  (MIS)"  for  the 
general  coordinate  geometry  used  in  this  analysis.  The 
variable  names  referenced  in  the  discussion  below  are  the 
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variable  names  used  in  the  non-MCNP  (user-supplied)  subrou¬ 
tines  or  functions. 

The  unmodified  version  of  MCNP  computes  distances  in  units  of 
distance,  centimeters.  However,  the  variable  density  atmos¬ 
pheric  model  computes  distances  in  units  of  mass  range, 
g/cm*  .  Therefore,  before  correction  for  a  variable  density 
atmosphere  can  be  made,  it  was  first  necessary  to  transform 
the  MCNP-computed  distance  to  collision,  PMF  (EDST  in  subrou¬ 
tine  EQDIST),  from  units  of  centimeters  to  an  equivalent  mass 
range  in  units  of  g/cm1  for  a  sea-level  homogeneous  atmos¬ 
phere.  This  equivalent  distance  was  then  corrected  for  a 
variable  density  atmosphere  using  the  MCNP  modifications. 

The  mass-integral  along  the  particle  path  in  a  homogeneous 
sea-level  atmosphere  is  calculated  using  the  following 
equation 


MIPDH-1 . 225x1  O^EDIST 


(2) 


where 


MIPDH  =  equivalent  mass  range  along  particle  path  in  a 
sea-level  homogeneous  atmosphere  [g/cm*  3 
EDIST  a  distance  to  collision  in  a  sea-level  homogeneous 
atmosphere  Ccm3 

This  calculation  is  performed  in  function  DMI . 

The  cumulative  mass-integral ,  MICA,  at  the  current  particle 
altitude,  ZO,  can  be  computed  using  Equation  (3)  below  once 
the  scale  height  region  has  been  determined.  This  calcu¬ 
lation  is  performed  in  function  MASI. 


where 


MICA  =  mass-integra 1  at  current  particle  altitude 
Cg/em*  3 

MI  Cs*.)  =  mass-integra  1  at  base  altitude  of  itr'  region 
[g/cm*  1 

p(z*>  =  density  at  base  altitude  of  i**'  region  Eg/c«*  3 
20  -  current  particle  altitude  Ecm3 
2*  =  base  altitude  of  i*"  region  Cc*3 
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S*.  =  mass-integral  scale  height  -factor  of  itr'  region 
Ccm3 

Zi  <  zO  <  Zi*i 

Though  the  collision  altitude,  ZC,  is  still  unknown,  the 
values  of  MIPDH,  MICA,  and  the  z-direction  cosine,  WO,  can  be 
used  to  determine  the  required  cumulative  mass-integral , 

MICP,  at  the  collision  altitude,  ZC.  MICP  is  computed  based 
on  the  known  mass  range  to  collision,  MIPDH.  This  calcu¬ 
lation  is  performed  by  calling  function  MIZ. 

MICP-MICA+MIPDH*  WO  (4) 


where 


MICP  =  required  mass-integral  at  collision  altitude 
Cg/cm2  3 

MICA  =  mass-integral  at  current  particle  altitude 
Cg/cm2  3 

MIPDH  =  mass-integral  along  particle  path  to  collision 
in  a  sea-level  homogeneous  atmosphere  Cg/cm2  3 

WO  =  z-direction  cosine 

At  this  point,  a  check  is  made  to  compare  the  computed  cumu¬ 
lative  mass-integral  at  the  collision  altitude,  MICP,  with 
the  mass-integral  at  infinity,  MI  INF,  defined  to  be 
1035.635131402448  g/cm2 ,  the  cumulative  mass-integral  at  the 
990  km  altitude  limit.  If  MICP  is  greater  than  MI  INF,  then 
Subroutine  EQDIST  is  exited  and  PMF  is  set  to  HUGE,  a  very 
large  number  used  by  MCNP,  effectively  eliminating  any  colli¬ 
sions  along  the  current  particle  track.  If  MICP  is  less  than 
or  equal  to  MI  INF,  then  the  collision  altitude,  ZC,  is  deter — 
mined  using  the  following  expression 

zc-Zj-sA4i+{"r^-MI’20)>| 

1  \  (p(Zt)Si)  J 


For  cases  where  the  z-direction  cosine,  WO,  becomes  very 
small  (nearly  co-altitude  relative  to  the  current  particle 
position),  the  change  in  density  is  very  small  along  the 
particle  path.  Therefore,  the  average  density  between  the 
current  altitude  and  the  collision  altitude  is  used.  In  this 
case,  the  difference  between  the  current  altitude  and  col¬ 
lision  altitude  is  <  10  meters.  The  air  densities  are  calcu¬ 
lated  by  calling  function  DNTY. 

Finally,  the  distance  to  collision,  PMF  (EDST  in  Subroutine 
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EQDIST),  corrected  for  a  variable  density  atmosphere,  can  be 
computed  by  calling  function  EQUIVDST.  Function  EQUIVDST 
accepts  as  inputs  the  collision  altitude,  ZC,  the  current 
particle  altitude,  ZO,  the  z-direction  cosine,  WO,  the  mass- 
integral  along  the  particle  direction  in  a  homogeneous  sea- 
level  atmosphere,  MIPDH,  and  the  densities  (calculated  only 
if  WO  is  small).  For  cases  where  WO  is  not  small,  the  cor — 
rected  distance  to  collision  is  given  by  Equation  (6)  below 

EDST-  (6) 


where 


EDST  =  distance  to  collision  corrected  for  a  variable 
density  atmosphere  [cm3 


For  cases  where  WO  is  small,  the  density  is  nearly  constant 
along  the  particle  path,  hence  the  corrected  distance  to 
collision  is  given  by 


msT- 


MIPDH 

[  (DENCA+DENCP)  /2] 


(7) 


where 


DENCA  =  air  density  at  current  particle  altitude  Cg/cm3] 
DENCP  =  air  density  at  collision  altitude  Cg/cm’J 

5.2.2  Mean  Free  Paths  To  Detector  or  Boundary 

Contributions  to  a  detector  are  made  at  every  source  or 
collision  point  by  creating  and  transporting  a  pseudoparticle 
directly  to  the  detector.  The  total  transmission  to  the 
detector  depends  on  several  factors:  1)  the  exponential 
attenuation  through  the  medium,  2)  a  probability  density 
function  for  scatter  toward  the  detector,  3)  spherical  diver¬ 
gence,  which  accounts  for  the  solid  angle  effect,  and  4)  the 
particle  weight. 

The  only  one  of  the  aforementioned  parameters  that  depends  on 
atmospheric  density  is  the  exponential  attenuation  term.  The 
attenuation  term  represents  the  total  attenuation  of  the 
radiation  by  the  atmosphere  along  the  particle  path  over  a 
given  number  of  mean  free  paths.  The  attenuation  along  the 
particle  path  to  a  detector  is  given  by  the  following 
relation 


e 


A-QXp(-AMFP) 


(8) 


where 


A  =  attenuation  along  the  particle  path 
AMFP  =  number  of  mean  -free  paths  to  the  detector  or 
nearest  intersecting  boundary 

At  this  point,  MCNP  computes  the  distance  to  the  detector, 

DD,  using  subroutine  DDDET .  MCNP  next  determines  the  minimum 
of  the  distances  to  the  detector  or  the  nearest  intersecting 
boundary.  MCNP  now  computes  the  macroscopic  cross  section  in 
the  cell,  PLE,  for  a  sea-level  homogeneous  atmosphere  using 
Equation  (1).  It  is  necessary  to  correct  the  macroscopic 
cross  section  in  the  cell,  PLE,  for  a  variable  density  atmos¬ 
phere  using  subroutine  EQDIST.  A  flag  is  set  prior  to 
calling  EQDIST  in  subroutine  TRANSM  to  indicate  the  origin  of 
the  calling  statement,  either  subroutine  HSTORY  or  subroutine 
TRANSM.  The  value  of  the  flag,  either  0  or  1 ,  will  determine 
the  specific  calculations  that  are  performed. 

Prior  to  calling  EQDIST,  the  reciprocal  of  the  macroscopic 
cross  section,  XMFP  =  1/PLE,  is  computed  to  give  the  mean 
free  path  in  a  sea-level  homogeneous  atmosphere  (recall  that 
MCNP  first  calculates  all  density-dependent  parameters  based 
on  sea-level  density  in  each  cell).  EQDIST  is  then  called 
from  subroutine  TRANSM  with  the  mean  free  path,  XMFP,  the 
current  altitude,  ZZZ,  the  z-direction  cosine,  WWW,  the 
calculation  flag,  NINP,  and  the  track  length  in  the  cell,  D, 
as  inputs. 

In  subroutine  EQDIST,  it  is  first  necessary  to  transform  the 
MCNP-computed  minimum  distance  to  the  detector  or  nearest 
intersecting  boundary,  D,  from  units  of  centimeters  to  an 
equivalent  mass  range  in  units  of  g/cm*  in  a  sea-level  homo¬ 
geneous  atmosphere.  This  equivalent  mass  range  is  then  used 
to  correct  the  macroscopic  cross  section,  PLE,  for  a  variable 
density  atmosphere. 

In  subroutine  EQDIST,  the  altitude  of  the  detector  or  nearest 
intersecting  bounding  surface  is  determined  using  the 
following  equation 


ZD-ZO+WO+DTD 


(9) 


where 


i 


ZP  *=  altitude  of  detector  or  intersecting  boundary  Ccm3 
ZO  =  current  particle  altitude  [cm3 
WO  =  2-direction  cosine 

DTD  =  distance  to  detector  or  intersecting  boundary  Ccm3 

rhe  mass -integral  along  the  particle  path  in  a  homogeneous 
atmosphere,  MIPDH,  is  next  calculated  using  Equation  (2). 
Function  DMI  performs  this  calculation. 

Ne.<t,  the  cumulative  mass-integral  at  the  current  particle 
altitude,  MICA,  and  the  cumulative  mass-integral  at  the 
detector  or  intersecting  boundary  altitude,  MIDA,  is  computed 
using  Equation  <3>.  This  calculation  is  performed  in 
function  MAS I . 

The  mass-integral  along  the  particle  path  in  a  variable 
density  atmosphere,  MIPD,  can  now  be  determined  using  the 
following  relation  (for  WO  not  too  small) 

Mipp.  no. 

wo 


where 


MIPD  =  mass-integral  ulong  the  particle  path  in  a 
variable  density  atmosphere  Cg/cm* 3 

For  cases  where  the  z-direction  cosine,  WO,  becomes  very 
small  (nearly  co-altitude  relative  to  the  current  particle 
position),  the  average  density  between  the  current  altitude 
and  the  collision  altitude  is  used.  In  this  case,  the  dif¬ 
ference  between  the  current  altitude  and  collision  altitude 
is  <  10  meters.  The  air  densities  are  calculated  by  calling 
function  DNTY.  For  small  values  o-.  the  z-di  reef  ion  cosine, 
WO.  MIPD  is  computed  using  the  fol 1  owing  relation 


MIPD-nTol 


(DENCA +DENCP) 


(U> 


where 

DTD  =  minimum  distance  to  the  detector  or  nearest  inter¬ 
secting  boundary  [cm3 

DENCA  =  air  density  at  current  particle  altitude  Cg/cm-l 
DENCP  =»  air  density  at  collision  altitude  Cg/cm33 
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Recall  that  the  minimum  of  the  distances  to  the  detector  and 
nearest  intersecting  boundary,  D,  was  first  transformed  into 
an  equivalent  mass  range,  MIPDH.  The  corrected  mass  range  in 
a  variable  density  atmosphere,  MIPD,  was  then  computed.  It 
is  now  possible  to  correct  the  mean  free  path  to  the  detector 
or  nearest  intersecting  boundary,  EDST,  for  a  variable 
density  atmosphere.  This  can  be  calculated  using 


ST-  (EDST*MIPDH) 
MIPD 


(12) 


where 


EDST  CLHS3  =  corrected  mean  free  path  CcmD 
EDST  CRHS3  =  uncorrected  mean  free  path  Ccml 
MIPDH  =  mass  range  along  particle  path  in  a  sea-level 
homogeneous  atmosphere  Cg/cm*  D 
MIPD  =  mass  integral  along  particle  path  in  a  variable 
density  atmosphere  Cg/cma ] 

Back  in  TRANSM,  the  corrected  macroscopic  cross  section  in 
the  cell,  PLE,  is  found  by  taking  the  reciprocal  of  the 
corrected  mean  free  path,  EDST  (defined  as  XMFP  in  MCNP), 
where  EDST  is  given  by  Equation  (12).  Finally,  the  cumu¬ 
lative  total  number  of  mean  free  paths  to  the  detector  or 
nearest  boundary  over  all  cells  traversed  by  the  particle 
trajectory,  is  found  from  the  following  equation 

AMFP-AMFP+PLE*D  <13> 


where 

AMFP  CLHS3  =  cumulative  total  number  of  mean  free  paths 
to  the  detector  or  nearest  boundary  along 
the  particle  path 

AMFP  CRHS3  =  number  of  mean  free  paths  along  the 
particle  path  in  the  current  cell 

AMFP  here  is  equal  to  PLE#D. 

6  User-Written  Subroutines  and  Functions 

This  section  describes  in  detail  the  implementation  of  the 
variable  density  model  developed  in  Section  5.  This  was 
accomplished  by  developing  various  user-written  subroutines 
and  functions  and  later  integrating  them  into  MCNP.  The 
variables  referenced  in  the  discussion  below  refer  to  the 
variable  names  used  in  the  user-written  subroutines  or 


•functions  as  opposed  to  the  variable  names  used  in  the  MCNP 
subroutines . 


6.1  Data  Blocks 

There  was  only  one  data  block  de-fined  as  part  of  the  MCNP 
modifications  and  is  described  below. 

6.1.1  Data  Block  ATINIT 

ATINIT  is  a  user — defined  data  block  which  initializes  the 
arrays  containing  parameters  which  define  the  variable 
density  atmospheric  model.  The  arrays  are  stared  in  a  user — 
defined  common  block  /ATCOM/.  Common  block  /ATCOM/  is 
defined  in  the  file  CM.inc. 

There  are  5  arrays  contained  in  data  block  ATINIT.  The 
arrays  SHFM  and  SHFD  contain  scale  height  factors  [cm]  for 
both  mass-integra  1  [g/cm*]  and  density  [g/cm3]  data,  respec¬ 
tively.  The  array  ZI  contains  the  base  altitude  [cm]  for 
each  scale  height  region.  The  array  DENI  contains  the 
reference  density  [g/cm3]  at  each  base  altitude.  The  array 
AMI  I  contains  the  cumulative  mass-integral  [g/cm2 ]  at  each 
base  altitude. 

6.2  Functions 

There  were  6  user — written  functions  incorporated  into  the 
MCNP  code.  The  function  of  each  of  these  is  explained  below. 

6.2.1  Function  DNTY 

This  function  calculates  a  density,  DNTY  [g/cm3],  given  an 
input  altitude,  Z  Ccm].  Function  DNTY  is  called  from  subrou¬ 
tine  EQDIST.  In  order  to  compute  a  density  based  on  the 
variable  density  atmospheric  model,  the  array  ZI  containing 
base  altitudes  for  each  scale  height  region,  the  array  DENI 
containing  densities  at  each  base  altitude,  and  the  array 
SHFD  containing  density  scale  height  factors  for  the  scale 
height  regions  are  required.  Function  DNTY  first  determines 
the  correct  scale  height  region  based  on  the  input  altitude, 
Z,  by  searching  the  base  altitude  array,  ZI.  The  array 
indexing  variable,  k,  is  used  to  index  the  DENI  and  SHFD 
arrays  to  obtain  the  corresponding  density  and  scale  height 
factor  for  that  scale  height  region.  Function  DNTY  then 
computes  the  density  at  Z  using  the  correct  atmospheric  model 
parameters. 

6.2.2  Function  MAS I 

This  function  calculates  a  mass-integral ,  MASI  Cg/cm2 ],  given 
an  input  altitude,  Z  [cm].  The  computed  mass-integral  repre- 
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sents  the  cumulative  mass-integral  -from  ground  level  up  to 
the  input  altitude,  Z.  Function  MASI  is  called  •from  subrou¬ 
tine  EQDIST.  In  order  to  compute  the  cumulative  mass- 
integral  based  on  the  variable  density  atmospheric  model , 
the  array  ZI  containing  base  altitudes  for  each  scale  height 
region,  the  array  DENI  containing  densities  at  each  base 
altitude,  the  array  SHFM  containing  density  scale  height 
factors  for  the  scale  height  regions,  and  the  array  AMII 
containing  the  cumulative  mass-integrals  for  the  base 
altitudes  are  required.  Function  MASI  first  determines  the 
correct  scale  height  region  based  on  the  input  altitude,  Z, 
by  searching  the  base  altitude  array,  ZI.  The  array  indexing 
variable,  k,  is  used  to  index  the  DENI,  SHFM,  and  AMII  arrays 
to  obtain  the  corresponding  density,  scale  height  factor,  and 
cumulative  mass-integral  for  that  scale  height  region. 
Function  MASI  then  computes  the  density  at  Z  using  the 
correct  atmospheric  model  parameters. 

6.2.3  Function  ZMAS 

This  function  calculates  an  altitude,  ZMAS  tern],  given  an 
input  mass-integral,  MI  Cg/cm*!.  The  computed  altitude 
represents  the  altitude  necessary  to  achieve  the  given  cumu¬ 
lative  mass-integral.  Function  ZMAS  is  called  from  subrou¬ 
tine  EQDIST.  In  order  to  compute  the  required  altitude  based 
on  the  variable  density  atmospheric  model,  the  array  71 
containing  base  altitudes  for  each  scale  height  region,  the 
array  DENI  containing  densities  at  each  base  altitude,  the 
array  SHFM  containing  mass-integral  scale  height  factors  for 
the  scale  height  regions,  and  the  array  AMII  containing  the 
cumulative  mass-integrals  for  the  base  altitudes  are 
required.  Function  ZMAS  first  determines  the  correct  scale 
height  region  based  on  the  input  mass-integral,  MI,  by 
searching  the  array,  AMII.  The  AMII  array  contains  the 
cumulative  mass-integral  for  the  base  altitudes  of  the  scale 
height  regions.  The  array  indexing  variable,  k,  is  used  to 
index  the  DENI,  SHFM,  and  ZI  arrays  to  obtain  the  corres¬ 
ponding  density,  scale  height  factor,  and  base  altitude  for 
that  scale  height  region.  Function  ZMAS  then  computes  the 
altitude  using  the  correct  atmospheric  model  parameters. 

6.2.4  Function  DMI 

This  function  calculates  a  mass-integral,  DMI  Cg/cm*3, 
between  two  points  in  a  sea-level  homogeneous  atmosphere. 

The  required  input  value  isi  1)  the  distance  to  collision  or 
mean  free  path  to  the  boundary  or  detector  [cm3 .  Function 
DMI  is  called  from  subroutine  EQDIST.  Function  DMI  computes 
the  mass  range  along  the  particle  path  of  known  range  in  a 
sea-level  homogeneous  atmosphere.  The  mass  density  at  sea- 
level  is  given  by  the  U.S.  Standard  Atmosphere,  1976  C33,  to 
be  1.225x10-*  g/cm*. 
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6.2.5  Function  MIZ 

This  function  calcuiates  a  mass-integral,  MIZ  Cg/cm23,  at  the 
new  particle  altitude  in  a  variable  density  atmosphere.  The 
new  particle  altitude  is  not  yet  explicitly  known.  The 
required  input  values  ares  1)  the  mass-integral  at  the 
current  particle  altitude,  MIP  Cg/cm23,  2)  the  mass  range 
along  the  particle  path  in  a  sea— level  homogeneous  atmos¬ 
phere,  DMI  Cg/cm*  3 ,  and  3)  the  2-direction  cosine,  W.  The 
z-direction  cosine  represents  the  cosine  of  the  angle  between 
the  polar  axis,  Z,  and  the  horizontal  .  is,  Y.  Function  MIZ 
is  called  from  subroutine  EQDIST.  Function  MIZ  computes  the 
mass-integral  required  at  the  new  particle  altitude  based  on 
previously  computed  mass-integral  values  and  the  current 
z-direction  cosine.  The  new  particle  altitude  is  computed 
subsequently. 

6.2.6  Function  EQUIVD5T 

This  function  calculates  an  equivalent  distance  to  collision, 
EQUIVDST  CcmD ,  along  the  particle  path  which  gives  the  same 
mass-integral  as  computed  in  a  sea- level  homogeneous  atmos¬ 
phere.  The  required  input  values  are:  i)  the  current 
particle  altitude,  Z  Ccm3,  2)  the  new  particle  altitude,  ZN 
Ccml,  3)  the  z-directian  cosine,  W,  4)  the  density  at  the 
current  particle  altitude,  DCA  Cg/cm3] ,  5)  the  density  at  the 
new  particle  altitude,  DCP  Cg/cm33,  and  6)  the  mass  range 
along  the  particle  path  in  a  sea-level  homogeneous  atmos¬ 
phere,  MIPD  Cg/cm*  3 .  Function  EQUIVDST  is  called  from  sub¬ 
routine  EQDIST.  If  the  difference  between  the  new  altitude 
and  the  current  altitude  is  less  than  1000  cm  (10  meters), 
then  W  is  very  small  and  the  density  changes  very  little 
along  the  particle  path.  In  this  case,  the  average  of  the 
densities  at  the  new  altitude  and  the  current  altitude  are 
used  in  computing  the  new  equivalent  distance  to  collision, 
EQUIVDST.  If  the  difference  between  the  new  altitude  and  the 
current  altitude  is  greater  or  equal  to  1000  cm  (10  meters), 
then  the  difference  in  the  altitudes  divided  by  the  direction 
cosine  is  used  to  compute  EQUIVDST. 

6.3  Subroutines 

There  was  only  one  subroutine  written  to  support  the  modifi¬ 
cations  built  into  MCNP.  This  subroutine  used  all  the 
functions  described  above.  This  subroutine  is  described 
below. 

6.3.1  Subroutine  EQDIST 

This  subroutine  is  used  to  calculate  one  of  two  transport  • 
parameters,  depending  on  the  calling  subroutine  within  MCNP. 
If  the  input  value  is  the  sampled  distance  to  collision,  PMF 
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[cm3,  -from  subroutine  HSTORY,  then  subroutine  EQDIST  returns 
an  equivalent  distance  to  collision  corrected  fa r  a  variable 
density  atmosphere,  EDST  [cm3.  If  the  input  value  is  the 
calculated  mean  -free  path,  XMFP  Ccm3,  to  the  detector  or 
nearest  intersecting  boundary  from  subroutine  TRANSM,  then 
subroutine  EQDIST  returns  an  equivalent  mear  free  path  cor¬ 
rected  for  a  variable  density  atmosphere.  If  subroutine 
EQDIST  is  called  from  subroutine  HSTORY,  then  the  input 
values  are:  1)  the  distance  to  collision,  EDST  [cm3 ,  2)  the 
current  particle  altitude,  ZO  [cm3,  3)  the  z-direction 
cosine,  WO,  4)  the  flag  NINP.  The  flag  NINP  has  either  the 
value  0  or  1 ,  indicating  which  calculations  are  to  be  per — 
formed.  NINP  has  the  value  0  if  the  calling  subroutine  is 
HSTORY,  and  hence  the  corrected  distance  to  collision  is 
required.  NINP  has  the  value  1  if  the  calling  subroutine  is 
TRANSM,  and  hence  the  corrected  mean  free  path  to  the 
detector  or  nearest  intersecting  boundary  is  required. 

Before  any  specific  calculations  are  performed,  subroutine 
EQDIST  first  defines  the  mass-integral  at  infinity,  MI  INF,  to 
be  1035.635131402448  g/cm* .  This  represents  the  cumulative 
mass-integral  at  the  upper  altitude  of  the  atmospheric  model, 
990  km.  Transport  calculations  are  not  performed  above  this 
altitude  limit. 

Next,  the  optimum  array  search  parameters  are  determined 
based  on  the  current  particle  location  and  direction.  Either 
the  upper  half,  the  lower  half  or  the  entire  array  is 
searched.  This  decreases  the  table  look-up  time  versus 
searching  the  entire  array  every  time.  The  array  search 
parameters  are  passed  as  input  values  to  the  six  usei — written 
functions  described  previously. 

At  this  point,  a  test  is  performed  to  check  the  value  of  the 
flag  NINP.  Recall  that  a  value  of  0  indicates  the  calling 
subroutine  is  HSTORY,  in  which  case  the  first  set  of  calcula¬ 
tions  are  performed.  These  calculations  are  primarily  a  set 
of  calling  statements  to  one  or  more  of  the  user — written 
functions.  Each  subsequent  line  of  code  is  well  documented 
as  to  its  function,  and  thus  will  not  be  repeated  here.  The 
theory  behind  each  calculation  is  described  in  Section  5  of 
this  addendum,  and  details  of  the  functions  performing  these 
calculations  are  described  above  under  the  heading 
"Functions*' .  It  is  important  to  note  that  before  the  new 
collision  altitude,  ZC  [cm3,  is  computed,  a  test  is  performed 
to  ensure  that  the  mass-integral  required  at  the  collision 
altitude,  MlCP,  is  within  both  the  upper  (MI  INF)  and  lower 
(0.0  g/cm* )  limits  of  the  atmospheric  model. 

If  this  test  is  successful,  then  calculations  continue.  The 
new  corrected  distance  to  collision,  EDST,  is  computed  and 
returned  to  subroutine  HSTORY.  In  subroutine  HSTORY,  PMF  is 


defined  to  be  EDST.  If  the  test  fails,  then  the  new 
corrected  distance  to  collision,  EDST  Ccm3,  is  returned  with 
the  value  -1.  In  subroutine  HSTORY,  EDST,  with  a  value  of 
-1,  is  defined  as  the  variable  PMF,  and  then  PMF  is  set  equal 
to  HUGE,  a  very  large  number  in  MCNP.  The  particle?  is  thus 
assumed  to  exit  the  problem,  effectively  eliminating  any 
further  particle  interactions  along  the  path. 

If  the  flag  has  the  value  1,  then  subroutine  EQDIST  is  called 
from  subroutine  TRANSM.  First,  the  detector  altitude  or 
altitude  of  the  nearest  intersecting  boundary  witn  the 
particle  path  is  computed.  The  given  input  values  of  the 
direction  cosine,  WO,  and  the  minimum  distance  to  either  the 
detector  or  the  nearest  intersecting  boundary,  DTD,  are  used. 
Each  subsequent  line  of  code  is  well  documented  as  to  its 
function,  and  thus  will  not  be  repeated  here.  The  theory 
behind  each  calculation  is  described  in  Section  5  of  this 
addendum,  and  details  of  the  functions  performing  these 
calculations  are  described  above  under  the  heading 
"Functions".  It  is  important  to  note  that  before  any  addi¬ 
tional  calculations  are  perf ormed ,  a  test  is  performed  to 
ensure  that  the  detector  or  boundary  crossing  altitude,  ZD, 
is  within  both  the  upper  (990  km)  and  lower  (0.0  km)  limits 
of  the  atmospheric  model.  If  this  test  is  successful,  then 
calculations  continue.  If  this  test  fails,  then  an  error 
message  is  printed  to  the  output  listing.  Calculations 
continue  until  a  new,  corrected  mean  free  path,  EDST  [cm3 ,  is 
computed.  EDST  is  based  on  the  given  input  mean  free  path, 
EDST  CcmD,  the  mass  range  along  the  path  in  a  sea-level 
homogeneous  atmosphere,  MIPDH  Cg/cm*3,  and  the  mass  range  in 
a  variaole  density  atmosphere,  MIPD  [g/cm* 3 .  EDST  is 
returned  to  subroutine  TRANSM  and  defined  to  be  XMFP,  the 
variable  used  in  subroutine  TRANSM. 

7  Modified  MCNP  Subroutines  and  Files 

This  section  describes  in  detail  the  auxiliary  modi f ications 
performed  on  various  MCNP  subroutines  and  files.  The  vari¬ 
ables  referenced  in  the  discussion  below  refer  to  the 
variable  names  used  in  the  MCNP  subroutines  or  files  as 
opposed  to  the  variable  names  used  in  the  user — written 
subroutines  or  functions. 

7.1  Include  Files 

There  were  two  include  files  that  were  modified,  and  these 
are  discussed  below. 

7.1.1  Ccmdeck  CM.inc 

This  include  file  serves  as  the  main  common  block  declaration 
file  for  all  overlays  used  in  MCNP.  The  following  modifica- 
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tions  to  CM.inc  were  mades 

1.  line  2:  the  -five  arrays  containing  parameters  for  the 
atmospheric  density  model  are  dimensioned  as  double  precision 
arrays . 

2.  line  37!  common  block  /ATCOM/  is  defined  and  contains  the 
necessary  arrays  that  store  the  various  parameters  defining 
the  variable  density  atmospheric  model. 

7.1.2  Comdeck  ZC.inc 

This  include  file  serves  as  an  auxiliary  comdeck  file  to 
comdeck  CM.inc  and  defines  processor-dependent  constants  and 
parameters  and  also  dimensions  general  constants  and  I/O  unit 
numbers.  The  following  modifications  to  ZC.inc  were  made: 

1.  line  6:  message  defined  to  inform  user  that  the  variable 
density  version  of  MCNP  is  being  used. 

2.  line  7:  error  statement  called  from  subroutine  EQDIST. 

3.  line  25:  two  variables,  NR1  and  NR2,  are  defined  and  used 
to  dimension  the  parameter  arrays  for  the  atmospheric  model. 
NR1  is  the  number  of  atmospheric  regions  and  NR2  is  the 
number  of  scale  height  factors. 

7.2  Subroutines 

There  were  four  subroutines  that  were  modified,  and  these  are 
discussed  below. 

7.2.1  Subroutine  IMCN 

This  subroutine  serves  as  the  main  overlay  subroutine,  and 
controls  the  input  and  sorting  of  the  problem  data.  The 
following  modifications  to  the  IMCN  overlay  were  made: 

1.  line  6  s  character  BLNK  defined  a^d  used  in  printing 
messages. 

2.  lines  112  thru  118!  message  printed  to  both  the  terminal 
and  output  listing  informing  the  user  that  the  variable 
density  version  of  MCNP  is  being  run. 

7.2.2  Subroutine  NEXTIT 

This  subroutine  controls  the  processing  of  input  items. 

1.  lines  34  thru  40:  density  in  all  cells  redefined  to  be 
I.225xl0-:5  g/em3,  the  mass  density  of  air  at  sea-level, 
regaroless  of  input  value. 
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7.2.3  Subroutine  HSTORY 


This  subroutine  performs  the  main  transport  calculations 
including  the  distance  to  collision,  PMF.  Starting  on  line 
61,  cross  sections  Cbarns/atoml  for  the  cell  are  determined 
via  subroutines  ACETOT  for  neutrons  and  PHOTOT  for  photons. 
On  line  75,  a  macroscopic  cross  section,  QPL  Ccm-1l,  is 
calculated.  In  the  variable  density  version  of  MCNP,  it  is 
not  recommended  that  the  exponential  transform  be  used  until 
a  full  study  is  made  of  the  effects  to  the  macroscopic  cross 
section.  The  mean  free  path,  GS  Ccml  is  next  computed  on 
line  81  by  taking  the  reciprocal  of  the  macroscopic  cross 
section,  QPL.  Up  to  this  point,  no  modifications  had  been 
made  to  subroutine  HSTORY.  The  following  modifications  were 
made  to  subroutine  HSTORY: 

1.  line  92:  random  sampling  of  the  logarithmic  distribution 
function  slightly  modified.  The  original  calculation  was  as 
fol lows 


qzzidg-RANGO 

(14) 

PMF— LOG  ( qzzidg)  +GS 

(IS) 

The  above  calculations  were  modified  as  follows 


qzzidg-RANGO 

(16) 

qzzidgl— LOG  (qzzidg) 

(17) 

PMF-qzzidgl *GS 

(IQ) 

where 

GS  =  mean  free  path  Ccml  based  on  the  macroscopic  cross 
section,  QPL 

PMF  a  distance  to  next  collision  Ccml 

The  purpose  of  this  modification  was  to  enable  the  recalcu¬ 
lation  of  the  mean  free  path,  SS,  and  the  macroscopic  cross 
section,  QPL.  These  values  then  serve  as  the  effective 
values  corrected  for  a  variable  density  atmosphere.  It  is 
possible  that  these  quantities  are  used  elsewhere  within 
MCNP,  affecting  quantities  not  directly  related  to  transport 
quantities.  A  full  investigation  of  this  was  not  made  due  to 
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time  constraints 


At  this  point,  transport  calculations  were  interrupted  in 
order  to  correct  the  above  quantities  for  a  variable  density 
atmosphere. 

1.  line  104s  the  flag  NINP  set  equal  to  0  indicating  the 
calling  subroutine  is  HSTORY. 

2.  line  105:  subroutine  EQDIST  called  to  calculate  the 
distance  to  next  collision,  PMF,  corrected  for  a  variable 
density  atmosphere. 

3.  line  108:  if  PMF  is  returned  with  the  value  -1,  then  PMF 
is  set  equal  to  HUGE,  a  very  large  number  used  by  MCNP.  This 
effectively  transports  the  particle  out  of  the  problem  with 
no  further  collisions. 

4.  lines  111  thru  113:  effective  values  of  the  mean  free 
path,  GS,  the  macroscopic  cross  section,  QPL,  and  the  macro¬ 
scopic  cross  section,  PLE,  are  recalculated.  Since  the 
exponential  transform  is  not  used  in  the  variable  density 
version  of  MCNP,  there  is  no  difference  between  QPL  and  PLE. 

Transport  calculations  continue  from  this  point  with  no 
further  modifications. 

7.2.4  Subroutine  TRANSM 

This  subroutine  calculates  the  transmission  to  the  detector 
or  nearest  intersecting  boundary  when  point  detectors  or  ring 
detectors  are  used.  Prior  to  calculating  the  transmission, 
subroutine  DDDET  is  called  to  calculate  the  distance,  DO,  to 
the  detector.  On  line  24,  subroutine  TRACK  is  called  to 
calculate  the  distance  to  the  nearest  intersecting  surface, 
DLS. 

Cross  sections  for  the  cell  are  determined  via  subroutines 
ACETQT  for  neutrons  and  PHQTOT  for  photons  on  1 ines  30  and 
31.  Once  the  microscopic  cross  section  Cbarns/atoml  in  the 
cell  has  been  determined,  it  is  multiplied  by  the  atom 
density  tatoms/bern-cm]  in  the  cell  to  obtain  the  macroscopic 
cross  section  tcm-13  on  line  32.  On  line  34,  the  track 
length  in  the  cell,  D,  is  determined  by  taking  the  minimum  of 
the  distance  to  the  nearest  bounding  surface,  DLS,  and  the 
distance  to  the  detector,  DD. 

Normally,  the  number  of  mean  free  paths  to  the  detector  or 
nearest  intersecting  boundary  is  computed  by  multiplying  the 
track  length  in  the  cell,  0,  with  the  macroscopic  cross 
section,  PLE.  At  this  point,  normal  calculations  are  inter¬ 
rupted  in  order  to  compute  a  new  macroscopic  cross  section  in 
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the  cell  corrected  for  a  variable  density  atmosphere.  The 
following  modifications  were  made  to  subroutine  TRANSM: 

1.  line  44:  the  flag  NINP  set  equal  to  1,  indicating  the 
calling  subroutine  is  TRANSM; 

2.  line  45:  the  mean  free  path,  XMFP  [cm]  is  calculated  by 
taking  the  reciprocal  of  the  macroscopic  cross  section,  PLE 
Ccm-1 ] ; 

3.  line  46:  subroutine  EQDIST  is  called  in  order  to  calculate 
a  corrected  mean  free  path,  XMFP; 

4.  line  48:  a  new  macroscopic  cross  section  is  determined  by 
taking  the  reciprocal  of  the  mean  free  path,  XMFP,  returned 
from  subroutine  EQDIST ; 

5.  line  51:  the  corrected  number  of  mean  free  paths  along  the 
track  length  to  the  detector  or  nearest  intersecting  boundary 
is  computed  using  the  corrected  macroscopic  cross  section, 
PLE. 

Calculations  continue  with  no  further  modifications. 

The  number  of  mean  free  paths  to  the  detector  or  nearest 
bounding  surface,  AMFP,  is  determined  using  the  values  of  the 
macroscopic  cross  section  in  the  cell,  PLE,  and  the  minimum 
distance  to  the  detector  or  nearest  bounding  surface,  D.  The 
locations  of  the  new  cell  and  particle  position  are  then 
updated . 

8  Miscellaneous  Information 

This  section  presents  some  miscellaneous  information  which 
might  be  useful  to  future  users  of  the  modified  MCNP  code  at 
AFIT.  Topics  discussed  include:  1)  the  nature  of  the  source 
spectra  used,  2)  the  location  of  the  MCNP  files,  and  3) 
comments  concerning  the  applicability  of  the  modified  MCNP 
code. 

8.1  Neutron-Photon  Source 

Thera  were  two  source  spectra  provided  with  the  SMAUG-1I  code 
C4],  The  spectra  included  37  neutron  energy  groups  ranging 
from  thermal  energies  to  a  maximum  of  14.2  MeV  and  21  photon 
groups  ranging  from  10“s  MeV  to  a  maximum  of  14  MeV.  The 
neutron  and  photon  spectrum  energies  and  corresponding  energy 
bin  probabilities  are  listed  in  Appendix  F  in  thesis  AFXT- 
/GNE/ENP/91M-6.  The  default  spectra  provided  with  the  SMAUG- 
II  code  were  used  in  the  MCNP  simulations  to  facilitate  the 
comparisons  between  SWAUG-II  and  MCNP.  Each  Monte  Carlo 
simulation  was  repeated  three  times  so  primary  neutrons, 
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primary  photons,  and  secondary  photons  could  be  generated. 
The  source  spectra  were  not  classified,  but  rather  generic, 
unclassified  spectra  provided  with  the  computer  code  SMAUG- 
II . 

8.2  Location  of  MCNP  files 

There  are  now  two  versions  of  the  MCNP  code  version  3B  at 
AFITs  1)  the  unmodified  version  of  MCNP  and  2)  the  modified 
version  of  MCNP  incorporating  a  variable  density  atmosphere. 
The  location  of  the  unmodified  MCNP  files  and  the  executable 
program  can  be  found  on  the  ELXSI  (GALAXY)  computer  at  AFIT, 
i.e.,  under  the  following  directory 

/ sre /men p/ Working 

The  location  of  the  modified  MCNP  files  and  the  executable 
program  can  be  found  under  the  following  directory 

/src/mcnp/Workingl 

The  location  of  the  cross  section  libraries  can  be  found 
under 


/src/mcnp/XCLib 

Copies  of  both  the  unmodified  and  modified  versions  of  MCNP 
wire  transferred  to  the  TRITON  SUN  system  at  AFIT,  in  LCDR 
Kirk  Mathews  office  (AFIT/ENP)  under  the  /dculp  and  /dmonti 
directories,  respectively .  Once  the  capability  to  run  MCNP 
on  TRITON  has  been  added,  either  MCNP  version  can  be  run. 
Copies  of  the  cross  section  libraries  were  also  transferred 
to  TRITON. 

The  file  transfer  sequence  to  copy  individual  files  from 
GALAXY  to  TRITON  are  given  below 

-  logon  to  userid  on  GALAXY 

-  go  to  directory  containing  applicable  files 

-  ftp  triton 

-  userid 

-  password 

-  cd  to  target  directory  on  triton 

-  mput  (multiple  files)  or  put  <1  file)  filename 

-  close 

-  quit 

If  a  large  number  of  files  needs  to  be  transferred  to  TRITON, 
use  the  following  sequence  of  commands 

-  logon  to  userid  on  GALAXY 

-  go  to  directory  containing  applicable  files 
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■  ar  r  archive  name  #.#  (archives  all  -files  in 
directory ) 

-  repeat  -ftp  sequence  above  to  trans-fer  the  archive  -file 
to  triton. 

Once  the  archive  file  has  been  transferred,  perform  the 
following  commands  to  retrieve  the  archived  files 

-  logon  to  userid  on  triton 

-  go  to  applicable  directory 

-  ar  x  archive  name  *.*  (extracts  all  archived  files) 
8.3  Comments 

8.3.1  Exponential  Transform 

Although  modifications  to  incorporate  a  variable  density 
atmosphere  to  MCNP  were  very  successful ,  further  work  is 
required  to  investigate  the  effects  of  the  exponential  trans¬ 
form  (not  used  thus  far)  on  the  macroscopic  cross  section. 

The  exponential  transform  can  be  invoked  in  the  HSTORY  sub¬ 
routine  from  the  input  file  (see  MCNP  user's  manual).  Since 
discrete  cells,  for  the  purpose  of  representing  the  variation 
of  air  density  with  altitude,  were  not  required,  all  trans¬ 
port  density-dependent  parameters  within  MCNP  were  first 
calculated  based  on  sea-level  air  density.  The  parameters 
were  then  corrected  for  a  variable  density  atmosphere.  The 
full  effects  of  invoking  the  exponential  transform,  which 
involves  path  length  stretching,  requires  investigation 
before  this  feature  can  be  used  with  confidence. 

8.3.2  Other  Density-Dependent  Features 

There  are  other  density-dependent  features  used  by  MCNP  which 
may  require  modif ication .  These  includes  1)  fission  heating, 
2)  dose-response  relationships ,  3)  cell  energy  deposition, 
and  4)  track  mean  free  path.  Fission  heating  (if  this  were 
somehow  invoked  in  the  atmosphere),  cell  energy  deposition, 
and  dose  calculations  all  involve  knowing  either  the  air 
density  and  the  density  of  some  other  material  (if  appli¬ 
cable)  within  the  atmosphere.  Modifications  would  be 
required  to  enable  the  user  to  use  the  variable  density 
version  of  MCNP  with  more  than  one  material.  The  track  mean 
free  path  is  printed  to  the  output  listing  upon  program 
completion.  It  is  computed  based  on  the  average  mean  free 
path  along  all  particle  tracks  within  a  cell.  Again,  fewer 
discrete  cells  are  used  in  the  modified  MCNP  code,  and  thus 
the  track  mean  free  path  within  a  cell  becomes  meaningless. 
This  will  not  affect  final  results;  the  track  mean  free  path 
is  used  only  as  diagnostic  information. 
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Appendix  A:  User-Written  Subroutines  and  Functions 


SUBROUTINE  EQDIST <EDST,  DTD,  ZO,  WO,  NINP) 


* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 


THIS  SUBROUTINE  CALCULATES  2  PARAMETERS,  DEPEND I NO  ON 
INPUT: 

1)  IF  THE  INPUT  VALUE  IS  THE  SAMPLED  DISTANCE  TO 
COLLISION  FROM  THE  ' HSTQRY .F '  SUBROUTINE,  THIS 
SUBROUTINE  RETURNS  AN  EQUIVALENT  DISTANCE  TO 
COLLISION  BASED  ON  A  VARIABLE  DENSITY  ATMOSPHERE. 

2)  IF  THE  INPUT  VALUE  IS  THE  CALCULATED  MEAN  FREE  PATH 
FROM  THE  ' TRANSM.F'  SUBROUTINE,  THIS  SUBROUTINE 
RETURNS  AN  EQUIVALENT  MEAN  FREE  PATH  BASED  ON  A 
VARIABLE  DENSITY  ATMOSPHERE. 

ALL  UNITS  ARE  IN  cgs,  CONSISTENT  WITH  MCNP. 


*=========  ■  - - == ,  ■  == — -  .  .  ■  == - 

* 

*  THE  FOLLOWING  IS  A  LIST  OF  ARRAYS  AND  VARIABLES  USED  IN 

*  THIS  SUBROUTINE  AND  EXTERNAL  TO  THE  FUNCTIONS. 

* 

*  — ■  ■  — 

* 

*  ZO 

t  wo 

*  EDST 

* 

*  DTD 

*  NINP 

t 

*  NR1 
$ 

t,  - - 

include  'CM.inc' 

EXTERNAL  MAS  I , DNTY , ZMAS , DMI ,MI 2 ,EQUIVDST 
DOUBLE  PRECISION  ZO, ZD , WO, MICA , DENCA , DENCP , ZC 
DOUBLE  PRECISION  MAS I , DNTY , ZMAS, DMI ,MI Z , EQUIVDST , DTD 
DOUBLE  PRECISION  EDST,MICP,MIDA,MI INF,MIPD,MIPDH 

C  DEFINE  MASS  INTEGRAL  AT  INFINITY  Cg/cnP  3  < ZO  >  990  KM) 

MI  INF  =  1035.635131402443 

C  DETERMINE  OPTIMUM  ARRAY  SEARCH  PARAMETERS 
L=NR1'2 

IFCWO.GE.O. .AND. ZO.LT.ZI(L) ) THEN 
1 1  =  1 


-  CURRENT  PARTICLE  ALTITUDE 

-  Z -DIRECT I ON  COSINE  RELATIVE  TO  POLAR  ANGLE 

-  EQUIVALENT  DISTANCE  TO  COLLISION  OR 

-  EQUIVALENT  MEAN  FREE  PATH  TO  BOUNDARY/DETECTOR 

-  DISTANCE  TO  BOUNDARY  OR  DETECTOR 

-  FLAB  INDICATING  WHICH  CALCULATIONS  ARE 
PERFORMED 

-  NUMBER  OF  SCALE  HEIGHT  REGIONS 
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non  on  on  on 


I2=NR1 

13=1 

ELSE I F  ( WO .  GE .  0 . .AND.ZO.GE. ZI (L) >THEN 
1 1  =L 
I2=NR1 
13=1 

ELSEIF(WO.LT.O. .AND.ZO.LT. Z I (L>  > THEN 
1 1  =L 
12=1 
1 3=- 1 

ELSE I F  <  WO . LT . 0 . .AND.ZO.GE. Z I <L) ) THEN 
I 1=NR1 
12=1 
1 3=- 1 

ENDIF 

I F ( N I NP . EQ . 0 ) THEN 


*  = - =====  ■  ■■  ■  .  == .  .  . — 

C  PERFORM  CALCULATIONS  BASED  ON  INPUTS  FROM  ' HSTORY . F ' 

C  SUBROUTINE. 

. ....  - - - - - . -  -  . 

C  CALCULATE  MASS  INTEGRAL  ALONG  PARTICLE  DIRECTION  IN  A 
C  HOMOGENEOUS  SEA-LEVEL  ATMOSPHERE. 

MIPDH  =  DMI(EDST) 

C  CALCULATE  MASS  INTEGRAL  AT  CURRENT  PARTICLE  ALTITUDE  IN 
C  A  VARIABLE  DENSITY  ATMOSPHERE. 

MICA  =  MMSitZO.il, 12, 13) 

CALCULATE  REQUIRED  MASS  INTEGRAL  AT  COLLISION  ALTITUDE  IN 
A  VARIABLE  DENSITY  ATMOSPHERE. 

MI CP  =  MIZ (MICA,  MIPDH,  WO) 

CHECK  TO  SEE  IF  COLLISION  ALTITUDE  IS  WITHIN  PROBLEM 
LIMITS. 

IF  (MICP.GE.O. .AND.MICP.LE.MI INF)  THEN 

CALCULATE  COLLISION  ALTITUDE  GIVEN  REQUIRED  MASS  INTEGRAL 
ALONG  PARTICLE  PATH. 

ZC  =  ZMASCMICP, II , 12, 13) 

CALCULATE  DENSITY  AT  CURRENT  PARTICLE  ALTITUDE  AND  AT 

COLLISION  ALTITUDE  FOR  CASES  WHERE  WO  IS  APPROX IMATELY 
EQUAL  TO  0. 

TF( (ZC-ZO) .LT.1000. >  THEN 
DENCA  =  DNTYtZO, II ,12,13) 

DENCP  =  DNTY (ZC, II, 12, 13) 

ENDIF 
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O  (J  **  (J  CJCJCJ  (JCJ 


C  CALCULATE  NEW  EQUIVALENT  DISTANCE  TO  COLLISION. 

EDST  =  EQUIVDST <  ZC,  ZO,  WO,  MIPDH,  DENCA,  DENCP ) 
ELSE 

EDST  =  -1. 

END  IF 

ELSEIF(NINP.EQ. 1 )THEN 


it . . . . .  .  . . 

PERFORM  CALCULATIONS  BASED  ON  INPUTS  FROM  'TRANSM.F' 
SUBROUTINE. 


CALCULATE  ALTITUDE  AT  DETECTOR  OR  BOUNDARY  CROSSING. 

ZD  =  ZO  +  WOtDTD 

IF ( ZD.LT. 0. . OR . ZD .GT . Z I (NR2) )THEN 
WRITE( IUO, ' ( 15X.A38) ' >MSG 
EDST=-1 
RETURN 
END  IF 

CALCULATE  MASS- INTEGRAL  ALONG  PARTICLE  DIRECTION  IN 
A  HOMOGENEOUS  ATMOSPHERE. 

SEA-LEVEL  ATMOSPHERE. 

MIPDH  =  DMI (DTD) 

CALCULATE  MASS  INTEGRAL  AT  CURRENT  PARTICLE  ALTITUDE  IN 
A  VARIABLE  DENSITY  ATMOSPHERE. 

MICA  =  MASHZO, II, 12,13) 

C  CALCULATE  MASS  INTEGRAL  AT  DETECTOR/BOUNDARY  ALTITUDE  IN 
C  A  VARIABLE  DENSITY  ATMOSPHERE. 

MIDA  =  MASKZD, 11,12, 13) 

C  CALCULATE  MASS  INTEGRAL  ALONG  PARTICLE  DIRECTION  IN 
C  A  VARIABLE  DENSITY  ATMOSPHERE. 

IF  (ABS(ZD-ZO) .GE. 1000)  THEN 
MIPD  =  (MI  DA— M I CA ) /WO 
ELSE 

DENCA  =  DNTY (ZO, II, 12, 13) 

DENCP  =  DNTY(ZD,U,  12,13) 

MIPD  *  DTD* ( DENCA+DENCP > /2 . 

END  IF 

C  CALCULATE  A  NEW  MEAN  FREE  PATH  BASED  ON  A  VARIABLE 
C  DENSITY  ATMOSPHERE. 

EDST  =  MIPDH  *  EDST  /  MIPD 


END  IF 

RETURN 

END 
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FUNCTION  EQUIVDST ( ZN,  Z,  W,  MIPD,  DCA,  DCP) 


*  ^ . . bb  — aa  - - 

* 

C  THIS  FUNCTION  CALCULATES  A  NEW  EQUIVALENT  DISTANCE  TO 
C  COLLISION  ALONG  THE  PARTICLE  DIRECTION  WHICH  GIVFS  THE 

C  SAME  MASS  INTEGRAL  AS  FOUND  IN  A  HOMOGENEOUS  SEA-LEVEL 

C  ATMOSPHERE . 

* 

*— - = . .  .  .  — - - - -  - 

......  .  -a -  -  - .  - 

* 

c  VARIABLES  USED  IN  THIS  FUNCTION: 

C 

C  ZN 

C 

C  Z 

C  W 

C 

C  DCA 

C  DCP 

C  MIPD 

C 

C  EQUIVDST 

C 

* 

». - -  •  . .  . . . 

DOUBLE  PRECISION  ZN, Z ,W, MIPD, DCA, DCP, EQUIVDST 
C 

C  IF  THE  DIFFERENCE  IN  ALTITUDES  BECOMES  SMALL  (  <  10m  ), 

C  USE  THE  AVERAGE  DENSITY  BETWEEN  THE  NEW  ALTITUDE  AND 
C  CURRENT  ALTITUDE. 

IF  (ABS(ZN-Z ) .GE. 1000)  THEN 
EQUIVDST  =  (ZN-ZJ/W 
ELSE 

EQUIVDST  =  MIPD  /  < 0. 5* < DCA+DCP ) > 

END  IF 
END 


FUNCTION  MAS I (Z, 11,12, 13) 

t, -  ■■  a; . -  _ _ gge . . . .  . . . . . . . . . . .  , .  

* 

GIVEN  AN  ALTITUDE,  FIND  THE  MASS  INTEGRAL. 

VARIABLES  USED  IN  THIS  FUNCTION: 

Z  -  GIVEN  PARTICLE  ALTITUDE 

ZI  -  ARRAY  CONTAINING  BASE  ALTITUDE  OF  REGIONS 
DENI  -  ARRAY  CONTAINING  DENSITY  AT  ZI 


-  NEW  ALTITUDE  IN  A  VARIABLE  DENSITY 
ATMOSPHERE 

-  CURRENT  PARTICLE  ALTITUDE 

-  Z-DIRECTION  COSINE  RELATIVE  TO  THE 
POLAR  ANGLE 

-  DENSITY  AT  CURRENT  PARTICLE  ALTITUDE 

-  DENSITY  AT  NEW  ALTITUDE 

-  MASS  INTEGRAL  ALONG  PARTICLE  DIRECTION 
(HOMOGENEOUS  SEA-LEVEL  ATMOSPHERE) 

-  EQUIVALENT  DISTANCE  TO  COLLISION  OR 
MEAN  FREE  PATH  TO  BOUNDARY/DETECTOR 
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AMI  I  -  ARRAY  CONTAINING  MASS  INTEGRAL  AT  ZI 
SHFM  -  ARRAY  CONTAINING  MASS  INTEGRAL  SCALE  HEIGHT 
FACTORS 

MAS  I  -  MASS  INTEGRAL  AT  GIVEN  ALTITUDE 

* 

*  . -  ■  ,  . - — - - - =============== 

include  'CM.inc' 

DOUBLE  PRECISION  MASI,Z 
DO  10  K=I 1 , 12, 13 

IF  (Z.GE.ZI  <K>  .AND.Z.LE.ZKK+l ) )  THEN 
MAS  I  = 

AMI  I  <K ) +DENI <K)*SHFM(K>*< 1-EXP(-(Z-ZI <K) )/SHFM(K) ) > 

END  IF 

10  CONTINUE 
END 


FUNCTION  ZMAS(MI , II , 12, 13) 


*  sa  ■  - . as ,  ■ - .  .  a  a, - -  —  — - ...  - —a——,  -  - —  — — 

* 

C  GIVEN  A  MASS  INTEGRAL,  FIND  THE  ALTITUDE. 

C 

C  VARIABLES  USED  IN  THIS  FUNCTION: 

C 

C  MI  GIVEN  MASS  INTEGRAL 

C  SHFM  -  ARRAY  CONTAINING  MASS  INTEGRAL  SCALE  HEIGHT 

C  FACTORS 

C  ZI  -  ARRAY  CONTAINING  BASE  ALTITUDE  OF  REGIONS 

C  AMII  -  ARRAY  CONTAINING  MASS  INTEGRAL  AT  ZI 

C  DENI  -  ARRAY  CONTAINING  DENSITY  AT  ZI 

C  ZMAS  -  ALTITUDE  CORRESPONDING  TO  GIVEN  MASS  INTEGRAL 

* 

*s==  a  eea  - - - sea - gaggcm - - 

include  'CM.inc’ 

DOUBLE  PRECISION  MI , ARG, ZMAS 
DO  30  K= 11,12, 13 

IF  (MI.GE.AMII(K).AND.MI.LE.AMII(K-U))  THEN 
ARG  =  LOG ( 1  +  ( AM 1 1  C  K  )  -  MI )  /  (DENI (K )  * 

SHFM ( K ) ) ) 

ZMAS  =  Z I (K )  -  SHFMIK)  *  ARG 
END  IF 

30  CONTINUE 
END 


n  n 


FUNCTION  DNTY( Z , 1 1 , 12, 13) 


*==  .  = - _= 

* 

C  GIVEN  AN  ALTITUDE,  FIND  THE  DENSITY. 

C 

C  VARIABLES  USED  IN  THIS  FUNCTION: 

C 

C  Z  -  GIVEN  PARTICLE  ALTITUDE 

C  ZI  -  ARRAY  CONTAINING  BASE  ALTITUDE  OF  REGIONS 

C  DENI  -  ARRAY  CONTAINING  DENSITY  AT  ZI 

C  SHFD  -  ARRAY  CONTAINING  DENSITY  SCALE  HEIGHT 

C  FACTORS 

C  DNTY  -  DENSITY  CORRESPONDING  TO  GIVEN  ALTITUDE 

* 

* - —  — ■  = . ===== — - - - =m ===== 

include  'CM.inc' 

DOUBLE  PRECISION  Z.DNTY 
DO  10  K=I 1 , 12 , 13 

IF  (Z.GE.ZI  <K>  .AND.Z.LE.ZHK+l  ) )  THEN 

DNTY  =  DENI <K>  *  EXP(-( Z-ZI (K) )  /  SHFD(K>) 

END  IF 

10  CONTINUE 
END 


FUNCTION  DMI (DIST) 


»x==B=r===as=  =  ■  =====  =«=  gam  . . . saces . . . r 

* 

C  CALCULATE  MASS  INTEGRAL  BETWEEN  TWO  POINTS  IN  A 

C  SEA-LEVEL  HOMOGENEOUS  ATMOSPHERE. 

C 

C  VARIABLES  USED  IN  THIS  FUNCTION: 

C 

C  DIST  -  DISTANCE  TO  COLLISION  OR  MEAN  FREE  PATH  TO 

C  BOUNDARY/DETECTOR 

C  DMI  -  MASS  INTEGRAL  ALONG  PARTICLE  DIRECTION 

t 

**■■■ sagas  asa  - - - - - am — gagas  ■  s  g  was as . ■-  a  c= = 

DOUBLE  PRECISION  DIST, DMI 

CALCULATE  MASS  INTEGRAL  ALONG  PARTICLE  PATH. 

DMI  »  1.22SE-3tDIST 

END 
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FUNCTION  MIZ (MIP,  DMI,  W) 


*  ==  =  ■  = - —  -  ..  ..  -  —-a-  - .  ■  - 

* 

C  CALCULATE  MASS  INTEGRAL  AT  NEW  ALTITUDE  IN  AN  EXPONENTIAL 
C  ATMOSPHERE. 

C 

C  VARIABLES  USED  IN  THIS  FUNCTIONS 

C 

C  MIP  -  MASS  INTEGRAL  AT  CURRENT  PARTICLE  ALTITUDE 

C  IN  A  VARIABLE  DENSITY  ATMOSPHERE 

C  DMI  -  MASS  INTEGRAL  ALONG  PARTICLE  DIRECTION  IN  A 

C  HOMOGENEOUS  SEA-LEVEL  ATMOSPHERE 

C  W  Z-DIRECTION  COSINE  RELATIVE  TO  THE  POLAR 

C  ANGLE 

C  MIZ  -  MASS  INTEGRAL  AT  NEW  ALTITUDE  IN  A 

C  VARIABLE  DENSITY  ATMOSPHERE 

* 

a -  -  - - -  ,  - 

DOUBLE  PRECISION  MIP , DMI , W , MI Z 

CALCULATE  MASS  INTEGRAL  AT  NEW  ALTITUDE. 

MIZ  =  W*DMI  +  MIP 
END 


BLOCK  DATA  AT  IN IT 


THIS  DATA  BLOCK  IS  USED  TO  SUPPORT  THE  VARIABLE 
DENSITY  ATMOSPHERIC  MODEL.  INITIALIZE  COMMON  /ATMCOM/ 
WITH  VALUES  OF  ATMOSPHERIC  PARAMETERS. 


include  ’CM.inc' 


* 

t 

* 


SHFM  -  ARRAY  CONTAINING  MASS  INTEGRAL  SCALE  HEIGHT 
FACTORS 

SHFD  -  ARRAY  CONTAINING  DENSITY  SCALE  HEIGHT  FACTORS 
ZI  -  ARRAY  CONTAINING  BASE  ALTITUDE  OF  VERTICAL 
REGIONS 

DENI  -  ARRAY  CONTAINING  DENSITY  AT  ZI 

AM 1 1  -  ARRAY  CONTAINING  MASS  INTEGRAL  AT  ZI 


DATA  SHFM/ 1 034 1 24 . 02237278 1 , 1 00998 1.9241 49543 , 

1  986693. 1726 19466, 964S03. 53272 1503B, 

2  941 155.9145785341,988092.312585042, 

3  692868.8443723,637442.9797604, 
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4 

5 

6 
7 

a 

9 

+ 

+ 

+ 

+ 

+ 

+ 


DATA 

1 

2 

3 

4 

5 

6 

7 

8 
9 
+ 


+ 

+ 


DATA 


i 


DATA 


DATA 


626 1 92 . 9959628 , 639736 . 906 17414, 

47557 . 78050254 , 7 1 0362 . 4275 1 1 288 1 , 

8382 1 1 . 26459594 1 5 , 776503 . 932 1 53373 , 
683277 .3511 724839 ,612608 . 545471 6644 , 
557270 . 926 1 50725 , 568864 . 8492265636 , 
619813.951 3853799 , 935 1 09 . 609 1771 976 , 
1267217 . 4 1 5038409 , 1 578526 . 59297959 , 
1870150. 85627943 , 22 1 9020 . 933057442 , 

27 1 4282 . 1 34795458 , 3257224 .61137614, 
3991584.358044663,4584848.61083246, 
5162659.761537237,5791446.057543592, 

550 1 266 . 555948576 , 6836368 . 56 1 027 1 77 , 
7834040 . 06307542 1,981 4372 . 605576273 , 

139 1 7805 . 4965649 , 1 8877525 . 6 1 508398/ 

SHFD/ 1 03039 1 . 726068488 , 1 006927 . 0550 1 6296 , 

983 1 53 . 5607496358 , 960534 .1411 264988 , 
937232 . 1 06253662 , 866352 . 1 969023093 , 
664086 . 7633898 1 56 , 637638 . 46529 , 

627630 . 265 119, 642604 . 7540969 , 

654589 . 3995756 , 7360 1 1 . 6057506 1 85 , 
834195.106975831,758287.3616637037, 
666098 . 1 02745308 , 592758 . 528 1 268486 , 
553227 . 469727463 1 , 5704 1 3 . 0947863342 , 

678 1 76 . 4068898259 , 997277 . 94 1 535646 , 
1324262.260093013,1632165.675462571 , 
1919412. 5389299 , 23 1 3392 . 496297635 , 

27974 1 2 . 470495623 , 340B28S . 286706647 , 
4112545. 8855796 1 6 , 4680457 . 5527 1 1 085 , 

53 11463. 853509 1 05 , 5890652 . 50294756 , 
6360111 . 158318114,6989325.564711994, 

8 1 62740 . 323944079 ,1041 9949 . 4345468 , 
14188168. 296 14109,1 9869924 . 24084625 / 

Z I /O . , 1 . D+5 , 2 . D+5 , 3 . 0+5 , 4 . D+5 , 5 . D+5 , 1 . D+6 , 

1 . 5D+6 , 2 . D+6 , 2 . 5D+6 , 3 . D+6 , 4 . D+6 , 5 . D+6 , 

6 . D+6 , 7 . D+6 , 0 . D+6 , 9 . D+6 , 1 . D+7 ,1.1 D+7 , 

1 .2D+7, 1 . 3D+7 , 1 . 4D+7 , 1 . 5D+7 , 1 .60+7,1 .8D+7, 
2 . 0D+7 , 2 . 4D+7 , 2 . 8D+7 , 3 . 2D+7 , 4 .0D+7 , 4 . BD+7 , 
5 . 6D+7 , 6 . 4D+7 , 7 . 2D+7 , 8 . 0D+7 , 8 . 8D+7 , 9 . 9D+7/ 
DEN 1 / 1 . 225D-3 ,1.111 7D-3 , 1 . 0066D- 3,9. 0925D-4 , 
8. 1935D-4 , 7 . 3643D-4 ,4.1 351 D-4 , 1 >  9476D-4 , 

B . 89 1 0-5 , 4 . 0084D-5 ,1.841 0-5 , 3 . 9957D-6 , 
l . 0269D-6 , 3 . 0968D-7 , 8 , 2829D-8 , 1 . 8458D-8 , 
3.4160-9,5. 604 D-i 0,9.7080-1 1,2.2220-11 , 

8. 1520-12,3. 83  ID- 12,2.0760-1 2, 1 .233D-12, 
5.1940-13,2.5410-13,7. 8580- 14,2.9710-14, 
1.2640-14,2. 803D- 15,7. 2080- 16,2. 0490- 1 6 , 
6. 523D-1 7,2.4480-17, 1 .1360-14,6.4640-18, 
3.7160-18/ 

AMI  I/O. , 1 16.7635,222.607167,318.334333, 

404.704533,482.4368,763.994467,91 1 .2723, 
978 . 75926666666669 , 1 009 . 379620 , 

1023 . 286286666667 , 1032 . 662946666667 , 
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1 034 . 006793333334 , 1 035 . 40648000000 1 , 
1 035 . 500609766667 , 1 035 . 624 1 07866667 , 
1 035 . 633205146667 , 1 035 . 634792366667 , 
1 035 . 635056 1 9600 1 , 1 035 . 635 1 04380667 , 
1 035 .63511 8027400 , 1 035 . 635 1 23665300 , 
1 035 . 635 1 26503 134,1 035 . 635 128111100, 
1 035 . 635 1 29736200 , 1 035 . 635 1 3047 1 237 , 
1 035 . 635 1 3 1 056504 ,1035. 635 1 3 1 2550 1 7 , 
1 035 . 635 1 3 1 334304 , 1 035 . 635 1 3 1 385704 , 
1 035 . 635 1 3 1 397859 , 1 035 . 635 1 3 1 400898 , 
1035.635131401864,1035.635131402191 , 
1035.635131 402325 , 1 035 . 635 1 3 1 402394 , 
1035.635131402448/ 


Appendix  B:  Modified  MCNP  Subroutines  and  Files 

C  INCLUDE  FILE  ZC.inc 

IMPLICIT  DOUBLE  PRECISION  (A-H.O-Z) 

C 

C  CODE  NAME  AND  VERSION  NUMBER. 

CHARACTER  KOD*8,VER*5 
PARAMETER  ( KOD=  *  MCNP ' , VER= ' 3B3 ' ) 

PARAMETER  ( MSG= ' DETECTOR  ALTITUDE  EXCEEDS  MODEL 
LIMITS' ) 

PARAMETER  <MSG1='***  NOTE:  USING  VARIABLE  DENSITY  MODEL 

***'  ) 

C 

C  PROCESSOR-DEPENDENT  NAMED  CONSTANTS. 

C  MDAS  IS  THE  INITIAL  SIZE  OF  DYNAMICALLY  ALLOCATED 

COMMON  C  /DAC/  ON  SYSTEMS  WHERE  MEMORY 

ADJUSTMENT  IS  NOT  AVAILABLE, 

C  SET  MDAS  LARGE  ENOUGH  FOR  YOUR  BIGGEST  PROBLEM. 

PARAMETER  ( MDAS=500000 ) 

PARAMETER  <NDP2=2 ,HUGE=1D37 ) 

PARAMETER  ( FTLS= . 2dO , DFT I NT= 1 00 . dO ) 

C 

C  ARRAY  DIMENSIONS.  I/O  UNIT  NUMBERS.  GENERAL 

CONSTANTS . 

PARAMETER 

( MAXE=SO , MAXF= 1 6 , MAX  V= 1 8 , MAX W=2 , MEMAX= 1 50 , M I NK-200 , 

1 

M I PT=2 , M J  SF=9 , MKFT =7 , MKTC=22 , MSEB=30 1 ,MXC=180,MXDT=20,MXDX=5, 
2 

MXLV= 1 0 , MXMTX= 1 00 , NBMX= 100 , NDEF= 1 4 , N0VR=5 , I U I ■ 1 , I UQ=2 , I UR=3 , 

3  IUX=4,IUD=7, IUB=8, IUP=9, IUS=10, IU1~1 1 , IU2=12, IUSW=13, 

4  IUSR=14, IUSC=15, IUC=16,IUT*17, IUZ=>18,  IUK=19,  JTTY=6, 

5 

ZER0=0 , dO , 0NE= 1 . dO , P I E=3 . 1 4 1 5926535B98d0 , F I VE 1 9=  < ZERO+5 . OdO ) * 
*19, 

6  AVGDN= . 59703 1 09d0 , GELEC= .511 OOBdO , GNEUT=939 . 58d0 , 

7  SL I TE=299. 7925dO , NR 1^36 , NR2S37 ) 

C 

C 

C - 

c 
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c 


INCLUDE  FILE  CM.inc 


include  'ZC.inc' 

DOUBLE  PRECISION 

SHFM ( NR )  , SHFD ( NR  > , Z I ( NR  > , DEN I ( NR ) , AM 1 1 ( NR ) 

C 

C  **********************  STATIC  COMMON 
************************* 

C 

C  FIXED  COMMON  —  CONSTANT  AFTER  THE  PROBLEM  IS 

INITIATED. 

COMMON  /FIXCOM/ 

ATWT(MEMAX) ,BCW<2,3> , DDG (M! PT , MXDT ) ,DXW(MIPT,3) , 

1 

DXX (MIPT , 5,MXDX ) ,  ECF (MIPT) ,EMCF(MIPT) ,EMX,ERGSAB(0:MAXE) , 

2 

ESPL(MIPT, 10)  ,FME( MEMAX  )  , FNW, RIM , RSSP , SNIT , SRV( 3 , MAXV ) , 

2 

TBLTMP(MAXE) ,TCO<MIPT) ,  THGF (0:50) ,WC1 ( MIPT ) ,WC2 ( MIPT ) ,WCSl (MI 
PT), 

3  WCS2 (MIPT) , WWG ( 7 ) , WWP (MI PT , 5 ) , 

3  ZFIXCM, 

4 

ICW, IDEFV(MAXV) , IDRC(MXDT) , IFFT, IGM, IGWW, IKZ , IMG , IMT , INK (MINK 

>, 

4 


IPLT,  I PTV ,  102,1  SSW ,  I  v’DD ( MAXF )  ,  I VD I S ( MAX V ) ,  I VORD ( MAXF )  , -JGM ( M I P 

T), 

5 

JTLX.JUNF, JXS(32,MAXE) , KFL , KNODS , KNRM , KPT (MIPT) ,KUFIL(2,6) , 

6 

KXS ( MAXE ) , LDR , LFCDG , LFCDJ , LME (MIPT, MEMAX ) , LMT ( MEMAX ) , LNP , 

6 

LOCDT ( 2 , MXDT ) , LVCDG , L  VCDJ , L  XS , MBNK , MCAL , MCT , MGEGBT (MIPT), MLAJ 
• 

7 

ML«JA,MODE,MRL,MSD,MSRK,MXA,MXAFM,NXAFS,MXE,MXF,MXFO,MXF2,MXF3 


* 

8 

MXJ ,MXT ,MXTR,NDET ( MIPT ) ,NDTT,NDX (MIPT ) ,NGWW(MIPT ) ,NISS,NJSR,N 
JSS, 

8 

NK  X  S  , NLEV , NL J  A , NNPOS , NOR  D , NP l , NP I KMT , NPN , NRCD , NRSS , NSPH , NSR , 

9  NSRCK , NTAL ,NTY (MAXE ) ,NVEC ,NWW( MIPT ) ,NXS( 16, MAXE) 

C 

C  OFFSETS  FOR  VIRTUAL  ARRAYS  IN  DYNAMICALLY  ALLOCATED 

STORAGE. 

COMMON  /FIXCOM/ 

LARA, LCMG (MIPT ) ,LDEN,LDXP,LEAA,LEWG(MIPT) ,LFIM, 

1 

LFMG,LF0R,LFRC,L6MG(MIPT) ,LGVL(MIPT) ,LGWT,LPMG(MIPT) ,LOAX,LRH 
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0,, 

2 

LSCF , LSMG <  M I PT ) , LSPF , LSQQ , LSSO , LTDS , LTHP , LTRF , LTTH , LVCL , LVEC , 

3 

LVOL , LWWE  (MIPT ) ,LWWF<MIPT> , LI  PA , LIPT , LI SS , LITD , LJAR , LJPT ,LJSC 

t 

4 

L J  SS , LJ  TF , L JUN , L J VC , LKCP  s  LKSD , LKST , LLAF , LL AT , LLCA , LLFC , LLFT , L 
LJA , 

5 

LLCT , LLST , LLSC , LMAT , LMFL , LMLL , LNCL , LNSF , LDDM , LDDN , LDEC , LDXC , L 
DXD, 

6 

LFLX(MIPT) , LFSO , LGWW  <  MIPT ) , LPAC , LPAN , LPCC , LPWB , LRKP , LTFC , LWNS 
* 

7 

LISE,LJFQ, LLAJ , LLCJ , LLSE , LNPW , LNSL , LNTB , LSCR , LDRC , LFDD , LGNR , L 
PIK, 

8  L I FL , L I GM , LPC2 , L JFL , L J  FT , LT AL , LBNK , LXSS , 

9  MFIXCM 
C 

C  THIS  IS  A  USER-DEFINED  COMMON  BLOCK  CONTAINING 

PARAMETER  ARRAYS 

C  FOR  THE  VARIABLE  DENSITY  ATMOSPHERIC  MODEL. 

COMMON  /ATMCOM/ 

SHFM ( NR 1 ) , SHFD ( NR 1 ) , Z I ( NR2 ) , DEN I ( NR2 ) , AM 1 1 ( NR2 ) 

C 

PARAMETER 

<NFIXCM=MIPT*MXDT+2*MEMAX+5*MIPT*MXDX+3*MAXV+2*MAXE+ 

1  25HCMIPT+71  ,LFIXCM= 

2 

3*MXDT+MINK+50*MAXE+< 1+MIPT ) *MEMAX+1 7*MIPT+2*MAXV+2*MAXF+159 ) 
DIMENSION  GFIXCM(NFIXCM) , JFI XCM ( LF I XCM ) 

EQU I VALENCE  <  ATWT , GF I XCM ) , ( I CW , JF I XCM ) 

C 

C  VARIABLE  COMMON  —  VARIABLE  BUT  REQUIRED  FOR  A 

CONTINUE  RUN. 

C  ARRAYS  THAT  ARE  BACKED  UP  WHEN  A  TRACK  IS  LOST. 

COMMON  /VARCOM/ 

FEBL  (2,  16  >  , PAX  <  MIPT  ,6,16),  RDUM(  50 )  ,  RLT ( 2  )  ,  SMIJL  (  3  )  t 

1  SUMK ( 3 ) , TMAV (MIPT  » 3  > , TWAC , TWSS , WSSI ( 7 ) , 

2  ZVARCM , 

3 

IDUM( 50  > , 1ST , IXAK , JRAD , KCSF , NBHWM, NBT ( MIPT ) ,NBY , NCT (MIPT ) , 

4 

NDRR ( MAXE ) , NETB  <  2  > , NPS , NQSW , NRSW , NSA f  NSS , NSS 1(8), NTSS , NWSB , NW 
SE, 

5  NWSG , NWSL , NWST , 

6  MVARCM 

PARAMETER  <NVARCM=99*MIPT+100, LVARCM=MAXE+2*MIPT+70) 
DIMENSION  GVARCM(NVARCM) , JVARCM ( LVARCM ) 
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EQU I  VALENCE  <  FEBL , GVARCM ) , < I DUM , JVARCM > 

C 

C  NOT-BACKED-UP  VARIABLE  COMMON 

COMMON  /NBVCOM/ 

CPK ,  CTS , DBCN  (  20 ) ,DDX (MIPT ,2 ,MXDX ) ,DMP , FIS , OSUM < 3 ) , 

1 

0SUM2(3,3> , PRN,RANI ,RANJ ,RIJK,RKK,RSUM(2) ,RSUM2<2,2> ,STLS, 

2  WGTS ( 2 ) ,  WTO, 

3  ZNBVCM , 

4 

KCT , KCY , KNOD , KTLS , KZKF , LOST  <  2 )  , NBOV , NERR , NFER , NLAJ , NLSE , NPC ( 2 
0), 

5 

NPD , NPNM , NPP , NPPM , NPSR , NQSS , NRN , NRRS , NSKK , NTC , NTC 1 , NWER , 

6  NWWS (2 , 99 ) , NZIP,NZIX,NZIY(MXDX , MI PT ) , 

7  MNBVCM 

PARAMETER  <NNBVCM=2*MIPT*MXDX+52 , LNBVCM=MIPT*MXDX+245 ) 
DIMENSION  GNBVCM (NNBVCM) , J NBVCM ( LNBVCM ) 

EQUIVALENCE  (CPK, GNBVCM) , (KCT,JNBVCM) 

C 

C  EPHEMERAL  COMMON  —  NOT  NEEDED  AFTER  THE  CURRENT 

RUN. 

COMMON  /EPHCOM/ 

ANG  <  3 ) , CPA , CTME , RANB , RANK , RANL , R ANS , SSB ( 10) , 

1 

TPP ( 20 ) , UDT (10,0: MXLV ) ,UDTS( 10*( 1+MXLV) ) , WNVP ( 4 ) , XHOM, YHOM , 

2  ZEPHCM, 

3 

I CHAN, ICS, IDMP, IFILE, ILN, ILN1 , IMTX(MXMTX ) , INDT , INFORM, IOVR, IT 
AL, 

4 

I TERM, ITFXS, ITOTNU, ITTY, IUOU , JCHAR , JFCN, JGF , JGXA ( 2 ) ,JGX0(2> , 

5 

JOVR ( NOVR ) , J VP , KONRUN , KPROD , LDQ , LFATL , LFLL , LGC (101), LSPEED , MI 
X, 

6  MNK , NBNK , NCH ( M I PT ) ,NDE , NKRP ,NST , 

7  MEPHCM 
PARAMETER 

( NEPHCM=66+20*MXLV , LEPHCM=MXMTX+N0VR+MIPT+1 37 ) 

DIMENSION  GEPHCM ( NEPHCM ) , JEPHCM ( LEPHCM ) 

EQUIVALENCE  (ANG, GEPHCM) , ( ICHAN, JEPHCM) 

C 

C  PBL  COMMON  —  PARTICLE  AND  COLLISION  DESCRIPTORS. 

C  IF  /PBLCQM/  IS  MODIFIED,  /PB9C0M/  MUST  BE  CHANGED  TO 

MATCH  IT. 

COMMON  /PBLCQM/ 

XXX ,YYY, ZZZ ,UUU, VVV,WWW,ERG, WGT,TME , VEL, I  CL , JSU, 

1 

IPT.NPA, I EX ,NODE , IDX ,NCP , KRN, JGP , DLS , DXL , DTC , F IML , F I  Ml ,FISMG, 

2  WTFASV , LEV , KKBNK , 1 1 1 , JJJ ,KKK, IAP. SPARE 1 , SPARE2 , SPARES , 

3  MPBLCM 

PARAMETER  (LPBLCM=NDP2*20+1 7) 


DIMENSION  J PBLCM ( LPBLCM ) , PBL ( 1 0 ) 

EQUIVALENCE  < XXX , JPBLCM,PBL> 

C 

COMMON  /TABLES/  EBL< 16) ,TALB<8,2> , JSF(MJSF) ,NVS(MAXV> 

C 

C  CHARACTER  COMMON  —  CHARACTER  VARIABLES  AND  ARRAYS. 

CHARACTER 

A I  D*80 , A I D1 *80 , AI DS*80 , CHCD* 1 0 , EXMS*80 , HBLN ( MAX V , 2  >  *3 , 

1 

HBLW ( M AX W ) *3 , HCS ( 2 ) * 7 , HFT  < MKFT ) *8 , HFU ( 2 ) *  11 , HMM ( MEMAX ) *  1 0 , 

2 

HMT  <  MXMTX ) * 10,HNP (MIPT ) *7 , HQVR*8 , HSD <  2  >*10,IBIN*8, IDTM* 19, 

3 

IDTMS* 19 , ILBL ( 8 ) *8 , KLIN *80 , KQDS*8 , KOVR ( NOVR ) *6 , KSF ( 29 ) *3, 

4 

LODDAT*B , L0DS*8 , MSUB < NDEF >*10, PROB I D* 1 9 , PROBS* 1 9 , RFQ (10) *57 , 

5  UFIL ( 3, 6 ) *  1 1 , VERS*5, XDATE < MAXE ) *8 , XLI ST ( MAXE ) *  10 , 

6  XSCRD ( MAXE ) * ( MXC ) 

COMMON  /CHARCM/ 

AID,AID1 , AIDS, CHCD , E  XMS , HBLN , HBLW , HCS , HFT , HFU , HMM , 

1 

HMT , HNP , HOVR , HSD , IBIN, IDTM, IDTMS, ILBL , KLIN ,KODS , KOVR , KSF , LODD 
AT, 

2  LODS , MSUB , PROB I D , PROBS , RFQ , UF I L , VERS , XDATE , XL I ST , X SCRD 

C  I SUB t  NAMES  OF  FILES 

CHARACTER*8 

I NP , OUTP , RUNTPE , SRCTP , XSD I R , P I X , WSSA , RSSA , COM , COMOUT , 

1  PLOTM , MCT AL , DUMN 1 , DUMN2 , I SUB ( NDEF ) 

COMMON  /CHARCM/ 

I NP , OUTP , RUNTPE , SRCTP , XSD I R , P I X , WSSA , RSSA , COM , 

1  COMOUT , PLOTM , MCTAL , DUMN 1 , DUMN2 
EQUIVALENCE  (ISUB.INP) 

C 

COMMON  /PB9C0M/ 

X  X  X9 , YYY9 ,1119, UUU9 , VVV9 , WWW 9 , ERG9 , WGT9 , TME9 , VEL9 , 

1 

I CL9 , J SU9 , IPT9 , NPA9 , IEX9, N0DE9 , IDX9, NCP9 , KRN9 , JGP9 , DLS9 , DXL9 , 

2 

DTC9 ,F IML9 ,F IM1 9 , FISMG9 , WTFAS9  , LEV9 ,i<KBNK9 ,1119,  JJJ9, KKK9 ,  LAP 
9, 

3  SP7,SP8,SP9, 

4  MPB9CM 

PARAMETER  < LPB9CM=LPBLCM> 

DIMENSION  JPB9CM ( LPB9CM ) 

EQUIVALENCE  ( XXX9 , JPB9CM ) 

C 

COMMON  /PB8C0M/ 

XXX8,YYYS , ZZZ8, UUU8 , VV V8 , WWW8 , ERG8 , WGT8 , TME8 , VEL8 , 

1  ICL8, JSU8, IPT8 
PARAMETER  ( LPB8CM=NDP2* 10+3 > 

DIMENSION  JPB8CM( LPB8CM > 
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EQUIVALENCE  < XXX8, JPB8CM) 

C 

COMMON  /BACKUP/  GVBU  <  NVARCM ) ,  J  VBU ( LVARCM ) 

C 

C  *****************  DYNAMICALLY  ALLOCATED  COMMON 
***************** 

C 

COMMON  /DAC/  DAS < MDAS/NDP2 ) 

C 

C  FIXED  DYNAMICALLY  ALLOCATED  COMMON. 

DIMENSION 

AAAFD(2) ,ARA( 1 ) , CM6( 1 ) , DEN< 1 > ,DXCP(MIPT, 1 ) ,EAA< 1 ) , 

1 

EWWG ( 1 ) , FIM(MIPT, 1 ) ,FMG< 1 > , FOR < MIPT, 1 ) ,FRC< 1 ) , GMG( 1 ) ,GVL( 1 ) , 

2 

GWT ( 1 ) ,PMG<1> ,QAX(MIPT,1),RH0<1 > , SCF < 1 > , SMG ( ' ) ,SPFC4,2> , 

3 

SQQ( 12,1) 4SSO( 1 ) , TDS ( 1 > ,TMP< 1 ) , TRF ( 17 , 0:0) , TTH( 1 ) ,VCL(3,7, 1 > , 

4  VEC (3,1), VOL <  1  >  , WWE ( 1 > , WWF ( 1 ) , 

5 

1 1 IFD< 1 ) , I PAN ( 1 > , IPTAL (8,5,1) ;I SS ( 1 )  , I TDS ( 1 ) , JASR ( 1 ) , JPTAL ( 8 , 
1), 

6 

JSCN( 1 ) , JSS( 1 ) , JTF(8, 1 ) , JUN( 1 > ,JVC< 1) ,KCP( 1 ) ,KSD<21 , 1 ) ,KST( 1 ) 

7  L AF (3,3) , L AT (2,1), LCA i 1 ) , LFCL (  1  )  , LFT ( MKFT ,  1  )  ,  L  J  A  (  1 )  , 

S 

LOCCT (MIPT , 1 ) , LOCST (MIPT , 1 ) ,LSC< 1 ) , MAT ( 1 ) ,MFL(3, 1 ! ,MLL(2, 1 ) , 

9  NCL(l) ,NSF< 1) 

EQUIVALENCE 

( DAS , AAAFD , ARA , CMG , DEN , DXCP . EAA , EWWG .  F I M ; FMR . FOR . FRC  t 
1 

GMG , GVL , GWT , PMG , QAX , RHO , SCF , SMG , SPF , SQQ , SSO , TDS , TMP , TRF , TTH , V 
CL, 

2 

VEC , VOL , WWE , WWF , 1 1 1 FD , I PAN , I PT AL , I SS , I TDS , J  ASR , J PT AL , J SCN , J  SS 
» 

3 

JTF, JUN,  JVC , KCP  , KSD , KST  , LAF , LAT  , LCA , LFCL  , LFT  , LJ A , LOCCT , LOCST  , 
LSC, 

4  MAT , MFL , MLL , NCL , NSF  > 

C 

C  VARIABLE  DYNAMICALLY  ALLOCATED  COMMON. 

DIMENSION 

AAAVD< 1 > ,DDM(2, 1 ) ,DDN(23, 1 ) ,DEC<2, 1 ) ,DXC(2, 1 > , 

1 

DXD(MIPT,23,MXDX) ,FLX( 1 ) ,FSO( 1 ) ,GWW(2,9, 1 ) ,PAC(MIPT , 10, 1 ) , 

2 

PAN(MIPT ,6,1 ) , PCC(3, 1 ) ,PWB(MIPT, 16, 1 ) , RKPL (5,1), TFC (3,20,1 ) , 

3  WNS (2,30), 

4 

II  IVD<1 )  ,  ISEF<2,  l )  , .JFQ(8,0: 1  >  ,LAJ(  1  )  ,LCAJ<  1  )  ,LSE(  1  >  ,NPSW<  1  >  , 
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5  NSL (10,1), NTBB ( 5 ,  i  ) 

EQUIVALENCE 

( DAS , AAAVD , DDM , DDN , DEC , DXC , DXD , FLX , FSO , GWW , PAC , PAN , 

1 

PCC , PWB , RKPL , TFC , WNS , 1 1 1 VD , I SEF , JFQ , LAJ , LCAJ , LSE, NPSW , NSL , NTB 
B) 

C  EPHEMERAL  DYNAMICALLY  ALLOCATED  COMMON. 

DIMENSION 

SCR ( 1 ) ,DRC< 16, 1 ) ,FDD<2, 1 > ,G£NR< 1 ) ,PIK< 1 ) , IFL< 1 > , 

1  IGMSAV (  1  )  ,  IPAC2 <  1  )  ,  JFL(  1  )  ,  JFK  1  ) 

EQUIVALENCE 

<DAS, SCR, DRC, FDD, GENR, PIK, IFL, IGMSAV, IPAC2, JFL, JFT) 

C 

C  TALLIES,  BANK,  AND  CROSS-SECTIONS  IN  DAC. 

DIMENSION  TAL ( 1 ) , I BNK ( 1  )  , XSS ( 1 ) 

EQUIVALENCE  ( DAS , TAL , I BNK , XSS) 

C 

C  CROSS  SECTIONS  ARE  REAL  ON  ALL  KINDS  OF  COMPUTERS. 

REAL  XSS 
REAL  YSS(l) 

EQUIVALENCE  (YSS.XSS) 

C 

C  DYNAMICALLY  ALLOCATED  COMMON  FOR  THE  IMCN  OVERLAY. 

DIMENSION  JTR (  1  )  ,AWTU  )  ,BBV(  1  )  ,PRB<  1  >  ,RTP(  1  )  ,SFB(  1  )  , 

1 

IPNK2.MKTC.0: 1 ) ,  JASW(  1 )  ,KAW(  1 )  ,KDUP(  1 )  ,KTR<  1 )  ,NLV(  1  >  ,NSLR< 10 

,D, 

2 

ARAS (2 , 1 ) , ATSA<  2, 1 > , RSCRN (2 , 1 ) , RSI NT (2,1), SCFQ (5,1), VOLS (2,1) 
* 

3  I INT  < 1 ) , I CRN ( 3 , 1 ) , LJAV ( 1 ) ,LJSV( 1 ) ,LSAT( 1 ) 

EQUIVALENCE 

( DAS , J  TR , AWT , BBV , PRB , RTP , SFB , I PNT , J  ASW , K AW , KDUP , KTR , 

1 

Nl ...  V ,  NSLR ,  ARAS ,  ATSA ,  RSCRN ,  RS I  NT ,  SCFQ ,  VOLS ,  1 1  NT ,  I  CRN  ,  LJAV.LJSV, 
2  LSAT ) 

C 

C  DYNAMICALLY  ALLOCATED  COMMON  FOR  THE  PLOT  OVERLAY. 

DIMENSION 

AMX  <4,4,1 ) , COE (6, 2, 1 ) ,CRS( 1 ) , JST(2, 1 > ,KCL( 102, 1 > ,KFM< 1 ) , 

1  LCL ( 1 ) , LSG  < 1 ) , NCS ( 1 ) ,PL8< 1 ) ,QMX<3,3,2, 1 ) ,ZST< 1 ) 
EQUIVALENCE 

( DAS , AMX , COE , CRS, JST ,KCL , KFM , LCL , LSG , NCS , PLB , QMX , ZST ) 

C 

C  DYNAMICALLY  ALLOCATED  COMMON  FOR  THE  MCPLOT  OVERLAY. 

DIMENSION 

AB1 ( 1 ) , AB2C 1 ) ,ERB( 1 ) .MCC< 1 ) ,ORD< 1 ) , XCC ( 1 ) ,YCC( 1 ) 

REAL  XRR { 1 ) ,YRR< 1 ) 

EQU I VALENCE  <  DAS , AB 1 , AB2 , ERB , MCC , ORD , XCC , XRR , YCC , YRR ) 

C 

c - 

C 
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SUBROUTINE  IMCN 

C  MAIN  CODE  OF  OVERLAY  IMCN. 

C  INITIATION  CODE  FOR  MONTE  CARLO  TRANSPORT, 

include  'CM.inc' 
include  'JC.inc' 

CHARACTER  BLNK*1 
C 

C  OPEN  SCRATCH  FILES  FOR  COLUMN  INPUT. 

HQVR= ' IMCN' 

OPEN ( IU1 ,STATUS=' SCRATCH' ) 

OPEN  < IU2, STATUS* ' SCRATCH ' ) 

C 

I F ( KQNRUN . NE . 0 ) GO  TO  190 
C 

***********************  INITIAL  RUN 
************************ 

* 

C  READ  THE  REST  OF  THE  INP  FILE  AND  SET  UP 

DYNAMICALLY 

C  ALLOCATED  STORAGE . 

DO  5  1=1 ,MINK*MNK 
5  INK ( I )=1 
CALL  PASS1 
HOVR= ' IMCN' 

I F  <  MODE . EQ . 0 ) KPT ( 1 >  =  1 

I F ( NSR . EQ . 6 . OR . ISSW. NE . 0 )CALL  SFILES 

I F  <  NJSR . EQ . 0 ) NJSR=NJSW 

IF ( NJSX . EQ . 0 ) NJSX=NJSR 

NDUP ( 3 ) =NDUP < 3 ) +NCPARF 

IF (NSR.EQ.71 . AND. NSRC.EQ.O) CALL  KSRCTP( 1 ) 

I F ( NSR . EQ . 7 1 ) MSRK=MAX ( 4500 , NSRCK+NSRCK/2 , NSRC/3 , MRL , MSRK ) 
ML J A=MLJ A+2*MX I T*NCOMP 
MLAJ=12*MXA+50 
IF ( MODE . NE . 2  >  MXT=MAX ( MXT , 1 ) 

IFtMODE.NE.l )NGWW(2-M0DE/2)=0 
CALL  SETDAS 

IF(LFLL.LT.LICC+4)CALL  CHGMEM (LFLL , L I CC+4 , ' IMCN  A 
C 

C  INITIALIZE  GENERAL  COMMON  NOT  YET  DONE. 

DO  30  1=1 ,MIPT 
DO  20  J=1 ,MXDT 
20  DDG  (  I ,  J  )  =HUGE 
DO  30  K= 1,2 
DO  30  J=1,MXDX 
30  DDX ( I , K | J ) =HUGE 
DBCN (10) =DFT I NT 
DMP=-240. 

ECF ( 2 ) = .00 1 
EMCF ( 2 ) =100* 

EMX=HUGE 

IKZ=5 
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DO  35  1=1 ,MAXF 
35  IVDD< I )=I 
KCY=1 
KTLS=1 
LOST ( 1 )=10 
LOST ( 2 ) = 10 
NPD=1000 
NTC=50 
NTC1 =50 
RKK= 1 . 

TOO ( 1 ) = . 00 1 *HUGE 
TOO ( 2 ) = . 00 1 *HUGE 
WC1 ( 1 ) =- . 5 
WC2 < 1 )=-  25 

I F ( MODE . NE . 0  >  WC 1 <2>=HUGE 
I F ( MODE . NE . 0 ) WC2 ( 2  >  =HUSE 
WGTS ( 1 )=HUGE 
WWP ( 1 , 1 ) =5 . 

WWP ( 2 , 1 ) =5 . 

WWP( 1 ,3)=5. 

WWP  ( 2 , 3 )  =5 . 

INITIALIZE  DYNAMICALLY  ALLOCATED  COMMON. 
DO  40  1=1 , (LICC+NDP2-1 1/NDP2 

40  AAAFD(I)=ZERO 

DO  41  I=LFCDG*NDP2+1,LFCDJ 

41  I I IFD ( I ) =0 

DO  43  I=LVCDG*NDP2+1,LVCDJ 
43  I I IVD ( I ) =0 

DO  45  I=LIFL+1 ,LAWT*NDP2 
45  IFL  < I )=0 

DO  47  I =L I PN+ 1 ,  L I CC 
47  JASW  < I ) =0 

DO  50  I=1,MXA 

IF ( NDX ( 1 ) . NE . 0 ) DXCP ( LDXP+ 1 , I )=1 . 

IF ( NDX  ( 2 )  .NE.0)DXCP(LDXP-*2,  I  >  =  1. 

50  IFtMODE.EQ. 1 )GWT(LGWT*I >=-l . 

IF ( NGWW ( 1 ) , NE . 0 ) EWWG ( LEWS ( 1 )+NGWW( 1)1=100. 

I F ( NGWW ( 2 ) . NE . 0 ) EWWG ( LT WG ( 2 ) ♦NGkW ( 2 ) >=100. 
LSC ( LLSC+1 ) =LSCF 
DO  60  I=0,NTAL 

60  IPNT(LIPN+1 , 1 , . ) =2+4+64+1 28 
TRF(LTRF+5,0)=1 . 

TRF  <  LTRF+9 , 0 ) = 1 . 

TRF ( LTRF+1 3 , 0 ) =1 . 

DO  65  J=1 ,M) TR 
DO  65  1=5,13 
65  TRF(LTRF+I , J)=HUGE 
DO  70  J=1 ,MXT 
DO  70  1=1, MXA 

70  TMP(LTMP+Ji-(  I-l  )*MXT)=253E-10 
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C  SORT  THE  TALLIES,  NEUTRONS  FIRST. 

DO  100  1=1 ,NTAL 
JFT (LJFT+I )=1 
K=  1 

DO  90  J=1,NTAL 

90  IF ( ABS(NTL ( J ) ) .LT . ABS < NTL CK > ) )K=J 
JPTAL (LJPT+l , I  )=MOD(NTL(K)  ,  1000) 

100  NTL <K >=5000 
C 

C  PRINT  MESSAGE  ON  VARIABLE  DENSITY  CODE  VERSION 

BLNK= '  ' 

WRITE (JTTY, *  >BLNK 

WRITE ( JTTY , ' (18x,A42) * >MSG1 

WRITE (JTTY , * ) BLNK 

WR I TE ( I UO , * ) BLNK 

WRITE ( IUO, ' ( 18x , A42 ) ' >MSG1 

WRITE < IUO,*) BLNK 

C  REREAD  AND  PRINT  THE  INP  MESSAGE  BLOCK  AND  TITLE 

LINE. 

REWIND  IUI 
DO  120  IP=1 ,3 
110  READ< IUI , * (A80) * >AID 
ILN=ILN+1 

WRITE (IUO, '  < 15 , 1H- , 7X , A80  > '  )ILN,AID 
IF( IP.EQ.2.AND.AID.NE. *  * )G0  TO  110 
120  IF(IP. EQ.l. AND. AIDC1*8>.NE. 'MESSAGE*' )G0  TO  130 
130  IF ( AID( 1 :5) .EQ. '  * . OR . A I D ( 6 : 72 ) . EQ . '  ' >G0  TO  140 

CALL  NXTSYM< AID, '  ',6,IT,IU> 

IF(KDATA(AID( 1 :5) ) . EQ . 2 . AND . KD AT A ( A I D ( IT : IU ) ) .EQ. 2 > 

1  CALL  ERPRNT (1,2, 0,0, 0,0, 0,0, 

2  ' 51HTHE  TITLE  CARD  LOOKS  SUSPICIOUSLY  LIKE  A  CELL 
CARD.  '  ) 

C 

C  REREAD  THE  REST  OF  THE  INP  FILE  AND  SET  UP  THE 

PROBLEM . 

140  CALL  RDPROB 

IF (MODE.EQ. 1 .AND. IFIP(2) .EQ . 0 . AND . NWW< 2 ) . EQ . 0 ) WRI TE ( IUO, ISO 
150  FORMAT ( /45H  PHOTON  IMPORTANCES  HAVE  BEEN  SET  EQUAL  TO 
l.  ) 

CALL  IGEOM 
CALL  ISOURC 
CALL  I  TALLY 
CALL  VOLUME 
HOVR=’ IMCN ' 

I F ( NSR . EQ . 40 ) CALL  WTCALC 
C 

C  WARN  OF  POSSIBLE  NEED  FOR  SUBROUTINE  SRCDX. 

DO  160  1=1 ,NTAL 

160 

IF ( JPTAL (LJPT+2, I ) .EQ . 5 . AND . JPTAL ( LJPT+3 , I ) . EQ . MAX ( MODE , 1 ) ) 
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1  GO  TO  170 

IF<NDX<MAX< MODE, 1 ) ) .EQ.O)GO  TO  180 
170  I F < NSR . EQ . 0 ) CALL  ERPRNT < 1,2, 0,0, 0,0, 0,0, 

1  ' 58HSUBR0UT I NE  SRCDX  IS  REQUIRED  IF  THE  SOURCE  IS 
ANISOTROPIC. * ) 

C 

C  WARN  OF  STRANGE  PHOTON  TIME  CUTOFF. 

180  I F < MODE . EQ . 1 . AND . TCO < 2 ) . NE.TCO( 1 ) >CALL 
ERPRNT  <1,2, 0,0, 0,0, 0,0, 

1  ' 55HPH0T0N  TIME  CUTOFF  IS  NOT  EQUAL  TO  NEUTRON  TIME 
CUTOFF . ' > 

C 

C  SET  UP  THE  WEIGHT  CUTOFFS. 

IF  < WC1 <  2 ) .EQ. HUGE. AND. MODE. EQ. 1 >WC1 <2>=WC1 < 1 ) 

1F(WC1 (2) . EQ.HUGE. AND. MODE. EQ. 2 >WC1 (2) =-.5 
I F  <  WC2  <  2  > . EQ.HUGE ) WC2  <  2 )  = . 5iWCl <  2 ) 

WCS1 < 1 ) =MAX  <  WC1 < 1 ) , -WC1 < 1 ) *SWTM ) 

WCS1 < 2 ) =MAX  < WC 1 ( 2 ) , -WC1 < 2) *SWTM> 

WCS2  < 1 ) =MAX  <  WC2 ( 1 > , -WC2  < 1 >  *SWTM ) 

WCS2  <  2 ) =MAX  <  WC2  <  2 ) , -WC2 ( 2 ) *SWTM  > 

IF  <NDE.NE.0)DBCN<2)=>NDE 
CALL  UFILES 
CLOSE < IU1 ) 

CLOSE  < IU2 ) 

RETURN 

C 

***********************  CONTINUE  RUN 
ft********:********  ****** 

* 

190  CALL  TPEFIL <  5 ) 

WRITE (IUO, ' < IX , A80) ' >AID 
I F < NSR . EQ . 6 . OR . ISSW.NE.OJCALL  SFILES 
DO  195  1=1 , MINK*MNK 
195  I NK  < I )  =  1 

I F < NSR . EQ . 7 1 > CALL  KSRCTP<3> 

IF<AID1 < 1 s Q ) .NE. 'CONTINUE' )G0  TO  240 
C 

C  REREAD  AND  PRINT  THE  INP  MESSAGE  BLOCK  AND  TITLE 

LINE. 

REWIND  IUI 
WRITE< IUO, ‘ < 1H  ) ‘ ) 

DO  210  IP=1 ,3 
200  READ< IUI , ‘ < AQO) ‘ )KLIN 
ILN=>ILN+1 

WRITE  <  IUO ,  *  <15,1  H- ,  7 X  ,  A80 )  *  ULN.KLIN 
IF < IP.EQ.2.AND.KLIN.NE. *  ')G0  TO  200 
210  IF< IP. EQ. 1 .AND. KLIN< 1 :9) .NE. 'MESSAGE: * )G0  TO  220 
C 

C  READ  THE  CONTINUE-RUN  DATA  FROM  THE  INP  FILE. 

220  CALL  RDPR08 
H0VR= ‘ IMCN* 

240  IF <NDE.NE.0)DBCN< 2 ) =ND£ 
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CALL  UFILES 
CLOSE < IU1 ) 

CLOSE < IU2) 

IF  * ( NFER . EQ . 0 . OR . LFATL . NE . 0 ) .AND. J0VR(4> .NE.O)CALL 
TPEFIL ( 6 ) 

RETURN 

END 


SUBROUTINE  NEXT IT 

C  PROCESS  THE  NEXT  INPUT  ITEM, 

include  'CM.inc' 

include  'JC.inc' 

CHARACTER  HT*75 
C 

C  CHECK  THE  ITEM  BEFORE  STORING  IT. 

NWC=NWC+1 

CALL  CHEKIT 

IF< I CS.LT.O) RETURN 

KS«INDEX ( ' ( >  ,HITM( 1 : 1 ) > 

C 

GO  TQ(  20,  60,  90,  95, 

10. 100.110.120.130. 140.150.161 . 165. 170.280, 

1 

180.190.210.216.216.210.218.220.280.280.280.280.280.230.280, 

2 

280 , 280 , 250 , 280 , 280 , 280 , 280 , 280 , 280 , 280 , 280 , 280 , 280 , 280 , 270 , 

3  280, 330, 340, 360, 370, 380, 390, 405, 410, 430) ICA 

GO 

TO  <  440 , 450 , 480 , 490 , 500 ,510, 520 , 530 , 540 , 550 , 560 , 570 , 580 , 590 , 60 
0, 

1  610,620,630,670,730,740,810,820,  10,930, 

2  1910, 1920, 1930) ICA-55 
10  RETURN 

C 

C  >>»>  CELL  DESCRIPTIONS 

C  M2CaPREV I OUS  SPECIAL  CHARACTER:  Q=NONE  1=(  2=) 

3=:  4=# 

C  M3C=FLAG  FOR  CELL  PARAMETERS. 

20  IFtHITM.EQ. ‘LIKE* .OR. L IKEF . NE . 0 )G0  TO  55 
I F ( KS . EQ . 0 . AND . K I TM . EG . 0  >  M3C= 1 
I F ( M3C . NE . 0 ) RETURN 
C 

C  STORE  THE  MATERIAL  NUMBER  AND  CELL  DENSITY. 

IF(NWC.EQ. 1 ) MAT (LMAT ♦MX A ) =1 ITM 

IF<NWC.EQ.2.AND.MAT(LMAT+MXA)  .NE.0)RH0(LRH0+t1XA)=RITM 
C 

*„■  - - - - - - - - - - - - - - -  .  =5  a  .  - -  -  egg . . . 

« 

C  SET  ALL  CELL  DENSITIES  TO  SEA-LEVEL  VALUE  Cg/cm=*3 

IF(NWC.EQ.2.AND.MAT(LMAT^MXA) .NE.O)  THEN 
IFtRITM.NE.O. )  THEN 
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on  nnnn  on  won  on 


RHO  <  L.RHO+MX  A )  =- 1 . 225E-3 
ELSE 

RHO ( LRHO+MX A ) =R I TM 
ENDIF 
END  IF 

* 

*  = 


I F ( NWC . EQ . 1 . OR . NWC . EQ . 2 . AND . MAT ( LMAT +MX  A ) . NE . 0 ) GO  TO  50 

PREPARE  TO  STORE  LOGICAL  OPERATOR  OR  SURFACE  NAME. 
NLJA=NLJA+1 
IF (KS.EQ.O)GO  TO  30 

STORE  LOGICAL  OPERATOR  AS  1000000+KS.  KS:  1  =  (  2=) 

=  :  4=# 

LCA <  LLCA+MX A ) =-ABS ( LCA ( LLCA+MX A ) ) 

LJ A ( LLJA+NLJA) = 1 OOOOOO+KS 
GO  TO  50 

STORE  THE  NAME  OF  A  SURFACE. 

30  LJA< LLJA+NLJA >=IITM 

50  M2C=KS 
RETURN 

55  CALL  LIKEBT ( 1  ) 

IF  (NWC.NE.  1  > RETURN 
DO  57  1=1, NDUP ( 1 ) 

57  IF  < KDUP ( LDUP+ 1 ) . EQ . I CN+ 1 00000  >  KDUP ( LDUP+ 1 ) =0 
RETURN 

>»>>  SURFACE  DESCRIPTIONS 

M1C=SURFACE  TYPE  INDEX. 

M2C=1  IF  SURFACE  TYPE  SYMBOL  IS  THE  SECOND  ITEM. 

60  I F ( K I TM . NE . 0 ) GO  TO  30 
M2C=NWC-1 
DO  70  M1C=1 , 29 

70  IF (KSF (M1C> . EQ. H I TM) RETURN 
80  IF  <  NWC . EQ . 1 ) JTR ( LJ TR+MX J )  =  1 1 TM 

I F ( NWC . GT . 1 ) SCF ( LSC (LLSC+MXJ ) +NWC-H2C- 1 >  =R I TM 
RETURN 

>>>>>  SPECIFICATION  OF  COORDINATE  TRANSFORMATIONS  FOR 
SURFACES  TR 

90  TRF (LTRF+1+NWC,MXTR) =RI TM 

I F ( NWC . L  T . 4 . OR . NWC . GT . 1 2  >  RETURN 

IF ( ICX . EQ. -1 ) TRF (LTRF+l+NWC ,MXTR ) =COS ( RI TMiPIE/ 1 60 . ) 

IF ( ABS<  TRF (LTRF+1 +NWC ,MXTR) -ANINT ( TRF ( LTRF+ 1 +NWC , MXTR ) ) ) .GT. 1 
E-10 ) 


1  TRF ( LTRF+ 1 , MXTR ) = I CN 
RETURN 
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>>>>>  VECTORS 
VECT 

95  N=(NWC+3>/4 

IF  ( NWC .  EQ .  4)KN-3 )  READ  (Hi  7M  (2:10)  , * (BN, 19) '  )  JVC<LJVC+N> 

I F  ( NWC .  NE  .  4*N-3 )  VEC  <  L  Vfc">NWC-4*N+3 , N  >  =R I TM 

RETURN 

>>>>>  CELL  IMPORTANCES 
IMP 

100  FIM ( LFIM+NQW , NWC ) =RITM 
IF( IFIP<2) .NE.O) RETURN 

I F ( MODE . EQ . 1 . AND . R I TM  r  GT . 0 . ) F IM ( LF IM+2 , NWC )  =  1 . 
IF(M0DE.EQ.2>FIM(LFIM+2,NWC)=FIM(LFIM+1 ,  NWC ) 

RETURN 

>>>>>  CELL  VOLUMES  FOR  TALLIES 
VOL 

1 10  I F  <  K I TM . NE . 0 ) VOL ( L VOL+NWC ) =R I TM 
I F  ( K I TM . EQ . 0 ) NOVOL= 1 
I F ( K I TM . EQ . 0 ) NWC=NWC- 1 
RETURN 

>>>>>  SURFACE  AREAS  FOR  TALLIES 
AREA 

120  ARA ( L ARA+NWC ) =R I TM 
RETURN 


>»>>  PHOTON  WEIGHT  LOWER  BOUNDS 
PWT 

130  GWT ( LGWT +NWC ) =R I TM 
RETURN 

>>>>>  EXPONENTIAL  TRANSFORM 
EXT 

140  1=1 

IF(HITM<ltl).EQ.  .OR.HITMU  s  1  )  .EQ.  ’ )I=2 
I F  <  H I TM ( I : I ) . NE . " S  * ) GO  TO  143 
J  =  I  +  1 
A=0 . 

GO  TO  149 

143  DO  i 45  J = I , N I TM 

145  IF ( INDEX  < ‘ VXYZ ' ,HZTM( J i J ) ) .NE.O) GO  TO  147 
147  HT=HITM( I i J-l > 

READCHT, *  <BN,E2i .0) * )A 
IF ( A . EQ . 0 . ) RETURN 
149  IF( J .EG.NITM+l )M=4 

IF(J.LT.NITM)READ(HITM(J  +  l i J ♦9 ) , ' (BN, 19)  •  )M 
IF  < J .LT.NITM)M=M+4 

I F ( J . EQ . N I TM ) M= I NDE  X (  ' XY2 ’  .HlTMCJsJ) ) 

QAX  ( LQAX+NQW ,  NWC )  =M+A 
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IF(HITMtlil) ,EQ. >  QA  X ( LQ AX +NQW , NWC ) =— M-A 
RETURN 

>>>>>  FORCED  COLLISIONS 
FCL 

150  FOR < LFOR+NQW ,NWC >  =RITM 
RETURN 

>>>>>  WEIGHT-WINDOW  LOWER  BOUNDS 
WWN 

161  WWF (LWWF ( NQW) + ( MAX (1,1 CN ) -1 ) *MXA+NWC )=RITM 
RETURN 

>>>>>  WEIGHT-WINDOW  ENERGIES 
WWE 

165  WWE ( LWWE ( NOW ) +NWC ) =R I TM 
RETURN 

>>>>>  WEIGHT-WINDOW  GAME  PARAMETERS 
WWP 

170  WWP ( NOW , NWC ) =R I TM 
RETURN 

>>>>>  DXTRAN  CELL  PROBABILITIES 
DXC 

180  DXCP(LDXP-*-NQW,NWC)=RITM 
RETURN 

>>>>>  CELLS  WHERE  FISSION  IS  TREATED  LIKE  CAPTURE 
NONU 

190  LFCL ( LLFC+NWC )  =  1 1 TM— 1 
RETURN 

>>>>>  SOURCE  DISTRIBUTIONS 
il.DS 

MALDISTRIBUTION  INDEX. 

M2C=CURRENT  LOCATION  IN  SPF . 

M3C=L0CAT I ON  OF  N  IN  KCP . 

M4C=NWC  OF  LAST  COLON. 

210  I F < KSD ( LKSD+20 ,  M  l  C ) .EG.O.OR,M2C.LT.KSD<LKSD*20,M1C>  >G0 
TO  213 

DO  212  IZ=M1C«-1  ,  MSD 
I=MSD«-M1C«-1-IZ 

KSD<LKSD«-13, I >=KSD<LKSD+13, I >+4 
M=MAX (KSDfLKSD+A, I ) , KSD  C LKSD+20 , 1 ) ) 

DO  212  J=1 ,M 
DO  212  K=l,4 

212  SPF<KSD(LKSD*13, 1 ) ♦K , M+ 1 - J ) =SPF  C KSD ( LK5D+ 1 3 , I ) ♦K,M-J ) 
MXXS“MXXS«-4 

213 

I F ( H I TM (1:1) . NE . ‘ D ' . OR . M4C . NE . 0 . AND . M4C . EQ . NWC- i . OR . N I TM . LT . 2 
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1  GO  TO  214 

IF (KDATA(HITM(2:NITM> > . EQ.O)GO  TO  214 

IF (KSD(LKSD+13,M1C) .EQ.0>KSD(LKSD+13,M1C)=MXXS 

M2C=M2C+1 

READ( HITM( 2 : 4  )  ,  '  <BN,E3.0> '  )SPF(KSD(LKSD-H3,M1C)  +  1  ,M2C> 
IF  ( KSD ( LKSD+6 ,  MIC ) . EQ . 0 ) SPF (KSD( LKSD+ 13 , MIC > +1 ,M2C>= 

1  SPF ( KSD ( LKSD+ 1 3 ,  M 1 C )  + 1 , M2C )  + 1 OOOOO 
RETURN 

214  IF (KITM.NE.O.OR.HITMt 1 t 1 ) .EQ. ' D' .OR.HITM.EQ. .OR. 

1  HITM.EQ. .OR.HITM.EQ. *)* >G0  TO  215 

NWC=NWC- 1 

IF ( INDEX ( ' LSFQT ' ,HITM< 1 : 1 ) ) .NE.0)KSD(LKSD+5,M1C>=1 

IF  (  INDEX  (  '  SQ '  , HI TM (  1 :  1  )  )  .  NE  .0)KSD  (LKSD+6,  Ml  0  =  1 

IF(HITM.EQ. 'Q' )KSD(LKSD+6 ,M1C>=1 

IF (HITM.EQ. ' T ’ ) KSD ( LKSD+9 , MIC > = 1 

IF (HITM.EQ. ' F ' )KSD(LKSD+1 1 ,M1C>=1 

IF (HITM.EQ. 'A' >KSD(LKSD+19,M1C>=1 

RETURN 

215  I F ( KSD ( LKSD+ 1 3 , M 1 C )  ..EQ .  0  >KSD(  LKSD+1 3 ,  MIC )  =MXXS 

IF ( ( M4C . EQ . 0 . OR . M4C . LT . NWC- 1 ) . AND.HITM.NE. * : ' .AND. 

1  HITM.NE. *C .AND.HITM.NE. ')' )GQ  TO  2157 

IF(HITM.NE. ' : * .AND.HITM.NE. * ( ' . OR .M4C. EQ .NWC-2 . AND.M4C .NE . 0 ) 
1  GO  TO  2155 
M3C=MKCP+2 
KCP ( LKCP+M3C- 1 )=-l 
KCP ( LKCP+M3C )  =  1 
MKCP=MKCP+3 

KCP  ( LKCP+MKCP )  =SPF  ( KSD  ( LKSD+ 1 3 , M 1 C>  + 1 ,  M2C ) 

SPF (KSD( LKSD+13 ,M1C ) +1  , M2C )  =-LKCP-M3C- 1 
2155  IF(HITM.EQ. ' : ' >M4C=NWC 
IF (HITM.EQ. ’ s ' ) RETURN 
IF  (HITM.EQ.  '  (  •  >M4C=NWC+3 
KCP ( LKCP+M3C ) =KCP ( LKCP+N3C  > + 1 
MKCP=MKCP+1 
KCP ( LKCP+MKCP ) « I  I TM 

IF (HITM.EQ. # ( • ) KCP (LKCP+MKCP >=1000001 
IF (HITM.EQ. ' ) ' ) KCP ( LKCP+MKCP  >  =  1000002 
IF(HITM( 1 i 1 > .NE. *D‘ ) RETURN 
READ ( H I TM ( 2 : 4 ) , * (BN, 13) ' > KCP (LKCP+MKCP) 

KCP ( LKCP+MKCP ) =KCP ( LKCP+MKCP ) ♦ 1 OOOOO 
RETURN 

2157  M2C=M2C+1 

SPF ( KSD ( LKSD+ 1 3 , M 1C ) + 1 , M2C ) «R I TM 
RETURN 

C  >>>>>  SOURCE  DISTRIBUTIONS 
SP,SB 

M 1 C=D I STR I  BUT  I ON  INDEX. 

M2C=FLAG  FOR  C. 

M3C=FLAG  FOR  FUNCTION. 

216  I = INDEX ( 4 PB ‘ , ICH( 2 : 2 ) )  ♦  1 
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I F  <  K I TM . NE . 0 ) GO  TO  217 
NWC=NWC-1 

IF  <  HITM. EQ . ' C ' >M2C=1 

IF (HITM.EQ. 'C' . AND.KSD(LKSD+19,M1C) . EQ.0)KSD(LKSD+19,M1C>=-1 
IF( HITM. EQ. 'V' ) KSD ( LKSD+10 , M 1C ) =1 
RETURN 

217  I F  ( NWC .  EQ .  1 .  AND .  1 1 TM ,  LT .  0 )  KSD  ( LKSD+2 , M 1 C )  =  1 1  TM 
I F ( NWC . EQ . 1 . AND . I ITM.LT.0)M3C=1 

I F ( M3C . NE . 0 ) SQQ ( LRQQ+NWC+3* I -6 , M 1 C  >  =R I TM 

IFCNWC.EQ.l .AND. IITM.EQ.-21 > SQQ (LSQQ+3* 1-4 , MIC > =12345 . 

I F ( M3C . NE . 0 ) RETURN 

IF ( KSD < LKSD+13 , MIC ) . EQ . 0 ) KSD (LKSD+13 , Ml C ) =MXXS 
SPF ( KSD ( LKSD+ 1 3 ,  M 1 C )  + 1  , NWC ) =R I TM 
RETURN 
C 

C  >>»>  SOURCE  COMMENT 
SC 

C  MALDISTRIBUTION  INDEX. 

218  HT=KLIN(6:80) 

DO  219  1=1,75,3 

219 

J3CN(MSSC+< I+2)/3)=ICHAR(HT( I : I > ) #65536+ I CHAR (HT ( 1+1 :  1  +  1 ) >*25 
6+ 

1  I CHAR ( HT  < I +2 : 1 +2  >  ) 

KSD ( LKSD+3 , M 1 C  >  =KSD  <  LKSD+3 , M 1 C  > +25 
MSSC=MSSC+25 
RETURN 
C 

C  >>>>>  SOURCE  DEFINITION 
SDEF 

C  M1C=NWC  OF  VARIABLE  NAME 

C  M2C=INDEX  OF  CURRENT  VARIABLE 

C  M3C=INDEX  OF  DEPENDED-ON  VARIABLE  OR  LOCATION  OF  N 

IN  KCP . 

C  M4C=NWC  OF  LAST  COlON. 

220  IF (HITM.EQ . *  =  *  >NWC=NWC-1 
IFIHITM.EQ. '  =  '  >  RETURN 

IF ( (HITM. EG. ' J ' .OR. HITM.EQ. '(' ) . AND . Ml C . EQ . NWC >MlC=MlC-2 

IF(M4C.NE.O. AND. MAC. EQ.NWC-2. AND. HITM. NE. ' : ' . AND .HITM. NE . ' < ' > 
1  M1C=NWC 

GO  TO ( 221 , 223 , 226 , 227 ) MIN ( 4 , NWC-M1 C+ 1 ) 

221  DO  222  M2C=1 ,MAXV 

222  I F ( H I TM . EQ . HBLN ( M2C , 1 > ) RETURN 

223  IF ( KITM .NE .0  >  GO  TO  227 
DO  224  M3C=1,MAXF 

224  I F ( H I TM . EQ . ' F ' / /HBLN  <  M3C , 1 ) >  GO  TO  225 
GO  TO  229 

225  I VDD ( M2C  >  =M3C 
RETURN 
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226  IF(M3C.NE.0. AND.HITM.NE. ' s ' .AND.HITM.NE. ' < ' )G0  TO  229 

227  IF(M4C.LT.NWC-1 .AND.HITM.NE. .AND.HITM.NE. '(' )60  TO 
228 

IF(HITM.EQ. ' : ' )M4C=NWC 

I F ( NWC . NE . M 1 C+2 ) GO  TO  2275 

M3C=MKCP+2 

KCP ( LKCP+M3C- 1 )=-l 

KCP ( LKCP+M3C  >  =  1 

MKCP=MKCP+3 

KCP  <  LKCP+MKCP ) =SRV  < 1 , M2C ) 

SRV ( 1 ,M2C)=-LKCP-M3C-1 

IF ( I VDIS ( 1 ) . NE . 0 ) KCP ( LKCP+MKCP ) = I VD I S  (  1  )  + 1 OOOOO 
I  VDIS ( 1 ) =0 

2275  IF<HITM.EQ. ' : ' ) RETURN 

IF<HITM.EQ. ' ( ' ) M4C=NWC+3 

KCP ( LKCP+M3C ) =KCP ( LKCP+M3C >  + 1 

MKCP=MKCP+1 

KCP  <  LKCP+MKCP  >  =  1 1 TM 

IF(HITM.EQ. ' ( ' ) KCP ( LKCP+MKCP ) =1000001 
IF(HITM.EQ. * ) '  ) KCP ( LKCP+MKCP  >  =  1000002 
IF(HITM(l:l) .NE. 'D' > RETURN 
READ(HITM(2:4 ) , * (BN, 13) ' ) KCP (LKCP+MKCP) 

KCP ( LKCP+MKCP ) =KCP ( LKCP+MKCP ) + 1 OOOOO 
RETURN 

228  SRV ( NWC-M 1C, M2C ) =R I TM 

I F ( NWC-M 1 C . LT . NVS ( M2C > ) RETURN 

229  M1C=NWC+1 
M3C=0 
RETURN 

C 

C  »>>>  TALLY  COMMENT 
FC 

230  HT=KLIN( 6:80) 

DO  240  1=1,75,3 

240  RTP(LRTP+IPL+(NWC-1 )*25+< I +2 ) /3 ) =ICHAR ( HT< I s I ) ) *65536+ 

1  I CHAR (HT( 1+1 s I+l ) )*256+ICHAR(HT< 1+2: 1+2) ) 

RETURN 

C 

C  >>>>>  ORDER  OF  TALLY  PRINTING 
FQ 

250  K= INDEX ( ' FDUSMCET ' ,HITM(1:1)) 

DO  260  1=1,7 

IF( JFQ(LJFO+I , ITAL) .EQ.K)JFQ(LJFQ+I , ITAL ) =JFQ ( LJFQ+I +1 , ITAL) 
260 

IF ( JFGKLJFO+I+l , ITAL) .EQ. JFOCLJFQ+I , ITAL) ) JFQ(L0FQ+I+1 , ITAL)= 
K 

JFQ ( LJFQ+9 , I TAL ) =-K 
RETURN 
C 

C  >>>>>  TALLY  FLUCTUATION  CHART  BINS 
TF 
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270  JTF(LJTF+NWC,ITAL>=IITM 
RETURN 
C 

C  >>>>>  OTHER  TALLY  CARDS 
PD, F ,FX ,FY ,FZ ,FT , E , T,C,FM, DE , 

C 

DF , EM , TM , CM , CF , SF , FS , SD , FU , DD 

C  M2C=NWC  OF  LAST  FT-CARD  PARAMETER  COUNT. 

C  SPECIAL  HANDLING  FOR  DXTRAN  DD  CARD. 

280  I F ( ICH.NE. ' DD ' . OR . I CN . GT . 2 ) GO  TO  300 
I =2-M0D ( NWC , 2  > 

J=(NWC+1 >/2 
IF ( I CN . EQ . 0 > GO  TO  290 
DDX  < I CN , I , J  > “R I TM 
RETURN 

290  I F  <  DDX ( 1 , 1 , J ) . EQ . HUGE ) DDX ( 1 , I , J ) =R 1 TM 
IF < DDX (2,1 , J ) . EQ.HUGE  >  DDX (2 , 1 , J  >=RITM 
C 

300  IF(KITM.EQ.O)GO  TO  310 
RTP ( LRTP+ 1 PL+NWC >  =R I TM 
RETURN 

310  I F ( I CH . EQ . ' FT ' ) GO  TO  323 
N=INDEX< 'NT' , HITM (1:1)  ) 

IF (N.EQ.O)GO  TO  320 

I  =  ( INDEX ( ‘FU  FS  FM  C  E  T  ' , ICH< 1 :3> )+5) /3 

I F ( I . EQ . 1 . AND . MOD  ( I  CN ,  1 0 ) .NE.5)I=0 

K=2**I 

L=MOD ( I PNT ( L I PN+ 1,1,1 T AL  > /K , 2 ) 

IF(L.NE.N-1 >IPNT<LIPN+1 , 1 , ITAL>=IPNT(LIPN+1 , 1 , ITAL)+(2*N-3)*K 
NWC=NWC~1 
RETURN 

320  I F ( KS . NE . 0  >  RTP  <  LRTP+ 1 PL+NWC )  =  1 OOOOOO+KS 
IF< ICH.NE. ' DE' .AND. ICH.NE. 'DF' ) RETURN 
IF (HITM. EQ, 'LIN' ' RTP ( LRTP+ I PL+NWC ) “- 1 . 

IF (HITM, EG. 'LOG' >NWC=NWC-1 
RETURN 

323  DO  325  I«1,MKF7 

325  IF(H1TM,EQ.HFT( I )  ) RTP  <  LRTP+ 1 PL+NWC  >  =  I 
I F ( M2C . NE , 0  >  RTP ( LRTP+ 1 PL+M2C ) =NWC-M2C- 1 
NWC"NWC+1 
M2C=NWC 
RETURN 
C 

C  >>>>>  DXTRAN  PARAMETERS 
DXT 

330  JcMOD( NWC-1 ,5)+l 
K=(NWC+4)/5 

IF(K.LE.MXDX)DXX(NQW,J ,K>=RITM 

IF(K.LE,tnXDX.AND.J  ,GE.4)DXX(NQW,  J  ,K)  =  (RITM*1 .00001  >**2 
IF ( J ,LE . 3 )DXW ( NQW , J ) =R I TM 
IF( J .NE.41RETURN 
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DXW(NQW, 1 >=0. 

DXW(NQW,2)=0. 

DXW(NQW,3>=0. 

RETURN 

C 

C  >>>>>  MATERIAL  SPECIFICATIONS 
M 

340  I F  <  MOD  <  NWC , 2 ) . EQ . 0  >  GO  TO  350 
MI X=MI X+l 
HT=HITM 

IF ( INDEX (HT,  '  .  '  )  .EQ.0)HT(NITM+1  :NITM+i  >  =  '  .  ' 

HMM<MIX ) (8- INDEX (HT, ' . * ) t 10>=HT 
I F ( HMM (MIX) (8:8) .NE. '  ' .AND. HMM (MI  X ) (9:9) .EQ. * 

' > HMM (MIX) ( 9  s  9  >  = ' 0 ' 

IF (HMM (MIX) (8! 10) .EG. '0  ' .OR. HMM (MIX) ( 8 : 10 ) . EQ . ' 00 

'  .OR. 

1  HMM (MIX) (8: 10) .EQ. '000' ) HMM (MIX) (8: 10)='  ' 

RETURN 

350  FME(MIX )=RITM 

I F ( R I TM . EQ . 0 . ) HMM ( M I X )  =  '  ' 

IF (RITM.EQ.O. )MIX=MIX-1 
RETURN 

C  >>>>>  NUCLIDES  FOR  DISCRETE  TREATMENT 
DRXS 

360  HT=HITM 

IF( INDEX (HT, ' . ' ) .EQ.0)HT(NITM+1 :NITM+1 ' 

HDR(NWC) ( 8- INDEX ( HT , * . * ) : 10)=HT 
RETURN 
C 

C  »>>>  TOTAL  OR  PROMPT  NUBAR 
TOTNU 

370  IT0TNU=2 
RETURN 
C 

C  >>>>>  ATOMIC  WEIGHTS 
AWTAB 

380  IF ( MOD < NWC , 2 ) .EQ. 1 ) KAW ( LKAW+ ( NWC+1 ) /2>  =  I ITM 
IF ( MOD  <  NWC  f  2 ) . EQ . 0  >  AWT ( LAWT+NWC/2  > =R I TM 
RETURN 
C 

C  >>>>>  CROSS-SECTION  DIRECTORY  INFORMATION 
XS 

C  M1C“LAST  CHARACTER  POSITION  USED  SO  FAR  IN 

XSCRD (NXSC ) . 

390  IF (NWC.NE. 1 )G0  TO  400 
NXSCs=NXSC+ 1 
XSCRD(NXSC )=' 

WR I TE( XSCRD ( NXSC ) (MXC-3:MXC) , ' ( 14) ' ) ICN 
XSCRD  <  NXSC  >  (8- INDEX  (HITM,  *  . '  )  ilO)<=HITM 
M1C=10 
RETURN 

400  XSCRD (NXSC) (M1C+2:M1C+NITM+1 )=HITM 
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MlOMlC+NITM+1 

RETURN 


>>>>>  VOID  CELLS 
VOID 

405  I0ID=0 

I=NAMCHG ( 1 , IITM) 

RH0(LRH0+I >=0. 

MAT (LMAT+I )=0 
RETURN 

>>>>>  PHYSICS  PARAMETERS 
PHYS 

410  I F ( NQW . EQ . 2 ) GO  TO  420 
IF (NWC.EQ. 1 )EMX=RITM 
I F ( NWC . EQ . 2 ) EMCF ( 1 ) =R I TM 
RETURN 

420  EMCF ( 2 ) =MAX ( R I TM , ZER0+ . 00 1 ) 

RETURN 

>>>>>  ENERGY  SPLITTING 
ESPLT 

430  ESPL ( NQW , NWC ) =R I TM 
RETURN 

>>>>>  THERMAL  TEMPERATURES 
TMP 

440  TMP ( LTMP+MAX  ( 1 , 1 CN  >  +  ( NWC- 1 ) *MXT  >  =R I TM 
RETURN 

>>>>>  THERMAL  TIMES 
THTME 

450  TTH ( LTTH+NWC ) =R I TM 
RETURN 

>>>>>  THERMAL  S(A,B>  DATA  SPECIFICATIONS 
MT 

480  INDT=INDT+1 

I MT  X ( I NDT )  =  I CN 
HT=HITM 

IF ( INDEX (HT, ' . •  ) . EQ . 0  >  HT ( NI TM+1 : NI TM+1 )  =  •  . * 

HMT ( INDT) (8-INDEX  (HT, ‘ .*  >$10)=HT 

I F ( HMT ( I NDT ) ( 8  s  8 ) . NE . ’  * . AND . HMT ( I NDT ) ( 9 : 9 ) . EQ . ‘  *) 
1  HMT( INDT) (9:9>='0‘ 

RETURN 

>>>>>  CUTOFFS 
CUT 

490  IF (NWC.EQ. 1 . AND , RI TM. NE . 0 . >TCO(NQW)=RITM 


IF (NWC.EQ. I .AND.RITM.NE.O. .AND.TC0<2) .EQ. . 001 tHUGE ) TCO( 2 ) =RIT 
M 

IF(NWC.EQ.2)ECF(NQW>«MAX ( ZERO ,RITM) 
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IF(NWC.EQ.3)WC1<NQW)=RITM 
I F ( NWC . EQ . 3  >  WC2 ( NQW )  =  . 5#  WC 1 ( NQW ) 

I F ( NWC . EQ . 4 ) WC2 ( NQW ) =S I GN ( R I TM , WC 1 ( NQW ) > 

I F ( NWC . EQ . 5 ) SWTM=R I TM 
I F ( NWC . EQ . 5 . AND . R I TM . EQ . 0 . >SWTM=-1 . 

RETURN 

>>>>>  SOURCE  PARTICLE  CUTOFF  NUMBER 
NPS 

500  NPP=I ITM 
RETURN 

>>>>>  COMPUTER  TIME  CUTOFF 

CTME 

510  CTME=RITM 
RETURN 

>>>>>  INTEGER  QUANTITIES  FOR  TEMPORARY  CODE  FEATURES 

IDUM 

520  I DUM ( NWC )  =  1 1 TM 
RETURN 

»>>>  REAL  QUANTITIES  ^OR  TEMPORARY  CODE  FEATURES 

RDUM 

530  RDUM ( NWC )=R ITM 
RETURN 

>>»>  PRINT  AND  DUMP  CONTROLS 

'RDMP 

540  I F ( NWC . EQ .  1 ) PRN=R ITM 
I F ( NWC . EQ . 2 ) DMP=R ITM 
I F ( NWC . EQ . 3  >  MCT = I ITM 
RETURN 

>>>>>  TERMINATION  AND  PRINT  CONTROL  FOR  LOST  PARTICLES 

LOST  550  LOST ( NWC >=I ITM 
RETURN 

>>>>>  DEBUGGING  CONTROLS 

DBCN 

560  DBCN < NWC )=R ITM 
RETURN 

>>>>>  SPECIFICATIONS  FOR  USER  FILES 
FILES 

570  J= (NWC+4 ) /5 
N=NWC-5* ( J-l ) 

IF(N.EQ. 1 )KUFIL(1 , J)=I ITM 
I F ( N . EQ . 2 ) UF I L ( 1 , J ) =H I TM 

IF ( N . EQ . 3 ) UF IL ( 2 , J ) =HSD( INDEX ( ' SD ‘  t HI TM (1:1))) 
IF(N.EQ.4)UFIL(3,J)=HFU( INDEX ( * FU* ,HITM< 1 1 1 >  >  > 

I F  <  N . EQ . 5 ) KUF IL(2,J>  =  IITM 
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RETURN 


>>>»  PRINT  CONTROL 
PRINT 

580  IF (MNK.EQ.l) RETURN 
IF ( I ITM.LT.01G0  TO  581 
INK< I ITM>-1 
RETURN 

581  I F ( MNK . NE . 0 ) GO  TO  585 
MNK=- 1 

DO  583  1=1, MINK 
583  INK( I )=1 
585  INK(-IITM>=0 
RETURN 

>>>>>  KCODE  SPECIFICATIONS 
CCODE 

590  I F  C  NWC . EQ . 4 ) KCT =R I TM 

I F ( NWC . EQ . 6 . AND . R I TM . NE . 0 . ) KNRM= 1 

I F ( NWC . EQ . 2 . AND . KONRUN . EQ . 0 . AND . R I TM . NE . 0 . )RKK=RITM 
I F ( NWC . EQ . 3 . AND . KONRUN . EQ . 0 ) I KZ=R I TM 
RETURN 

>>>>>  LOCATIONS  OF  KCODE  SOURCE  POINTS 
KSRC 

600  FSO ( LFSO+NWC+2* ( (NWC-l >/3> )=RITM 
RETURN 

>>>>>  WEIGHT-WINDOW  GENERATOR  PARAMETERS 
WWG 

610  WWG ( NWC  >  =R I TM 
RETURN 

>>>>>  ENERGY  BINS  FOR  WEIGHT-WINDOW  GENERATOR 
WWGE 

620  EWWG ( LEWG ( NQW ) +NWC ) =R I TM 
RETURN 

>>>>>  SURFACE  SOURCE  WRITE  INFORMATION 
SSW 

M1C=NWC  OF  VARIABLE  NAME. 

M2C= INDEX  OF  CURRENT  VARIABLE. 

630  I F ( NWC . GT . N J  SS ) GO  TO  640 
M1C=NJSS«-1 
JSS ( LJSS+NWC  >  =  1 1 TM 
RETURN 

SEARCH  FOR  KEYWORDS  AFTER  ALL  SURFACE  NUMBERS  ARE 

READ . 

640  IF(HITM.EQ. '=• )NWC=NWC-1 
IF (HITM.HQ. •=• > RETURN 
I F ( K I TM . NE . 0 ) GO  TO  660 
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M1C=NWC 

DO  650  M2C=1,MAXW 
650  I F ( H I TM . EQ . HBLW ( M2C ) ) RETURN 
660  I F ( M2C . EQ . 1 ) NSPH= 1 1 TM 
IF(M2C.EQ.2> IPTY=I ITM 
RETURN 

>>>>>  SURFACE  SOURCE  READ  INFORMATION 
SSR 

M1C=NWC  OF  VARIABLE  NAME. 

M2C=INDEX  OF  CURRENT  VARIABLE. 

670  IF<HITM.EQ.  '  ■=' )NWC=NWC-1 
IF(HITM.EQ. ) RETURN 
I F ( NWC . GT . M 1 C ) GO  TO  690 
DO  680  M2C=1,MAXV 
680  IF (HITM .EQ .HBLN(M2C , 2 ) ) RETURN 
690  I F  <  K I TM . EQ . 0 ) GO  TO  720 
I F ( M2C . GT . 3  >  GO  TO  710 
I F ( M2C . GT . 2 ) GO  TO  700 
J  ASR ( L  J  AR+NWC-M 1 C  >  =  1 1 TM 
I F ( NWC-M 1C.LT.NJSR) RETURN 
GO  TO  720 

700  I SS  <  L I SS+NWC-M 1 C ) = 1 1 TM 

I F ( NWC-M 1 C . LT . N J  S  X  >  RETURN 
GO  TO  720 

710  SRV ( NWC-M 1C, M2C ) =R I TM 

I F ( NWC-M 1 C . LT . NVS <  M2C > ) RETURN 
720  M1C=NWC+1 
RETURN 

»>>>  UNIVERSE  DESIGNATORS 
U 

730  J  UN ( L J  UN+NWC )  =  1 1 TM 
RETURN 

»>»  CELL  TRANSFORMATIONS 
TRCL 

M1C=CELL  INDEX 

M2C=NWC  OF  LAST  CELL  SEEN,  NEGATIVE  IF  LEFT  PARENS. 
740  GO  TO (750, 770, 780) 1+INDEXC * ( ) * ,HITM( 1 : 1 ) ) 

750  IF(M2C.LT.0)G0  TO  760 
M 1 C=M 1 C+NWC— M2C 
KTR<LKTR+M1C)=IITH 
M2C=NWC 
RETURN 

760  M=NWC+M2C+1 

TRF  <LTRF+M,MXTR)=RITM 
I F ( M . LT . 5 . OR . M . GT . 13) RETURN 

IF ( ICX.EQ.-l ) TRF( LTRF+M , MXTR ) =COS ( R 1 TMtPIE/ 1 80 . ) 

IF ( ABS( TRF (LTRF+M ,MXTR ) -AN  I NT ( TRF ( LTRF+M , MXTR ) ) ) .GT. IE-10) 

1  TRF(LTRF  +  l , MXTR)  =>ABS(  TRF (LTRF+1  , MXTR)  ) 
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RETURN 

770  M 1 C=M 1 C+NWC-M2C 
MXTR=MXTR+1 

KTR  ( LKTR+M  10  =  1  OOO+MXTR 
TRF (LTRF+1 ,MXTR)=-1000-MXTR 
M2C=-NWC 
RETURN 

780  I F  <  NWC . NE . -M2C+2 ) GD  TO  790 

KTR ( LKTR+M 1 C ) =TRF ( LTRF+2 , MXTR ) 

TRF ( LTRF+2 , MXTR  >  =HUGE 

MXTR=MXTR-1 

GO  TO  800 

790  CALL  TRFMAT ( MXTR  > 

800  M2C=NWC 
RETURN 

>>>>>  LATTICE  TYPE 
LAT 

810  I F  < 1 1 TM . EQ . 0 ) RETURN 
LAT ( LLAT + 1 , NWC )  =  1 1 TM 
NLAT=NLAT+1 
LAT ( LLAT+2 f  NWC  >  =NLAT 
RETURN 

>>>>>  CELL-FILLING  UNIVERSES,  WITH  TRANSFORMATIONS 
FILL 

M1C=CELL  INDEX 

M2C ( POS I T I VE ) =NWC  OF  LAST  ITEM  SEEN,  NOT  IN  PARENS 
M2C ( NEGAT I VE ) =-NWC  OF  LEFT  PARENS,  WHEN  IN  PARENS 
M3C=0  WHEN  NOT  IN  LATTICE  FILL 
M3C ( NEGAT I VE ) =-NWC  OF  12  IN  11:12 
M3C ( POS I T I VE ) =N  OF  LAF ( LLAF+M , N ) 

M4C=MAX IMUM  VALUE  OF  N 

820  GO  TO ( 830 , 880 , 890 , 900) 1+INDEX (':()' , HI TM( 1:1)) 

830  I F ( M2C . LT . 0 ) GO  TO  870 
I F ( M3C . EQ  *  0 ) GO  TO  860 
I F ( M3C . GT . 0 ) GO  TO  840 
I =NWC+M3C 

IF (MODI  1,3) . NE . 0 )LAF (LLAF+MLAF+2+I /3 , 1 )  =  I ITM 
IF ( MODI  I , 3) . EQ . 0 ) LAF (LLAF+MLAF+l+I/3, 2 )=I ITM- 
1  LAF (LLAF+MLAF+l+I/3, 1 >+l 
IF ( I .NE. 6) RETURN 
M3C=2 
M2C=NWC 

M4C=2+LAF ( LLAF+MLAF+ 1,2) #LAF ( LLAF+MLAF+2 , 2 ) * 

1  LAF ( LL AF  +ML AF +3,2) 

RETURN 

840  M3C=M3C+NWC-M2C 

I p ( M3C . GT . M4C ) GO  TO  850 
LAF < LLAF+MLAF+ 1 , M3C )  =  1 1 TM 
M2C=NWC 
RETURN 
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850  MLAF=MLAF+M4C*3 
M3C=0 

860  M1C=M1C+NWC-M2C 

MFL(LMFL+1 ,M1C>=IITM 

M2C=NWC 

RETURN 

870  M=NWC+M2C+1 

TRF ( LTRF+M , MXTR ) =R I TM 
IF(M.LT.5.0R.M.GT. 13)RETURN 

IF ( ICX.EQ.-l ) TRF  <  LTRF+M, MXTR ) =C0S ( RITMfcPIE/ 180 . ) 

IF(ABS(TRF ( LTRF+M , MXTR > -AN I NT (TRF ( LTRF+M , MXTR ) > ) .GT. IE-10) 
1  TRF ( LTRF+1 , MXTR ) =ABS ( TRF ( LTRF+1 , MXTR > ) 

RETURN 

880  IF (M3C.NE.0) RETURN 
M3C=-NWC-1 

LAF ( LLAF+MLAF+ 1,1) =MFL ( LMFL+ 1 ,  M 1 C ) 

MFL ( LMFL+ 1 ,  M 1 C  >  =-LL AF-ML AF 
RETURN 

890  MXTR=MXTR+1 

IF ( M3C . EQ . 0 ) MFL  (  LMFL+3 ,M1C )  =1000+MXTR 

I F ( M3C . NE . 0 ) LAF  <  LLAF+MLAF+3 , M3C )  =  1 OOO+MXTR 

TRF ( LTRF+ 1 , MXTR  >  =- 1 OOO-MXTR 

M2C=-NWC 

RETURN 

900  I F ( NWC . NE . — M2C+2 ) GO  TO  910 

IF ( M3C .EQ .0 ) MFL ( LMFL+3, M 1C  >=TRF (LTRF+2 , MXTR ) 

I F ( M3C . NE . 0 ) LAF  < LLAF+MLAF+3 , M3C ) =TRF ( LTRF+2 , MXTR ) 

TRF ( LTRF+2 , MXTR ) =HUGE 

MXTR=MXTR-1 

GO  TO  920 

910  CALL  TRFMAT ( MXTR  > 

920  M2C*NWC 
RETURN 

>>>>>  PHOTON-PRODUCTION  BIAS 
PIKMT 

930  NPI KMT  =NP IKMT  + 1 

PIK ( LP IK+NPIKMT ) =RI  TM 
RETURN 

>>>>>  FIRST  SPARE  CARD  TYPE 
2A 

1910  CONTINUE 
RETURN 

>>>>>  SECOND  SPARE  CARD  TYPE 
ZB 

1920  CONTINUE 
RETURN 

>>»>  THIRD  SPARE  CARD  TYPE 
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zc 

1930  CONTINUE 
RETURN 
C 

END 


SUBROUTINE  HSTORY 

C  RUN  THE  COMPLETE  HISTORY  OF  A  SOURCE  PARTICLE, 

include  'CM.inc' 
include  'RC.inc' 

C 

C  DEBUG  FEATURES:  SET  UP  EVENT  LOG.  PRINT  DEBUG 

LINE. 

KRFLG=0 

I F ( NPS+ 1 . GE . DBCN ( 3 ) . AND . NPS+ 1 . LE . DBCN ( 4 )  ) KRFLG= 1 
I F ( DBCN ( 2 ) .EQ.O. ) GO  TO  20 

IF (M0D(NPS+1 , I NT (DBCN (2) > ) .EQ.O)WRITE( IUO, 101NPS+1 ,NCT< 1 >+NCT 
(2), 

1  NRN.RIJK 

10  FORMAT ( 5H  DBCN, 319, 4X ,F16.0,TL1 , 1H  > 

C 

C  SAVE  VARIABLE  COMMON  FOR  POSSIBLE  LOST  PARTICLE 

RERUN . 

20  DO  25  1=1 ,NVARCM 
25  GVBU  Cl) =GVARCM  < I ) 

DO  30  1=1 ,LVARCM 
30  JVBU( I )=JVARCM( I ) 

C 

C  START  A  PARTICLE  FROM  THE  SOURCE. 

40  CALL  STAR TP 

I F ( NTER . NE . 0 ) GO  TO  310 
I F ( KDB . NE . O ) GO  TO  410 
IF(NST.NE.O)  then 
RETURN 
end  i-f 
C 

C  TERMINATE  THE  PARTICLE  IF  ITS  ENERGY  IS  BELOW 

CUTOFF . 

C  MOST  BANKED  PARTICLES  COME  BACK  HERE. 

60  I F ( ERG . LT . ECF ( I PT ) . EQ V . MCAL . LT . 2 ) GO  TO  270 
C 

C  CALCULATE  THE  DISTANCE  TO  THE  CELL  BOUNDARY,  DLS. 

I F ( LCA ( LLCA* I CL  >  » LT . 0  >  CALL  CHKCEL Cl  CL , 3 , J  ) 

70  I F ( WGT . L£ . 0 . >  CALL  EXP  I  RE ( 1 , * HSTORY ‘ , 

l  ‘THE  WEIGHT  OF  THE  CURRENT  PARTICLE  IS  ZERO  OR  LESS.') 
CALL  TRACK (I CL) 

I F ( KDB . NE . 0 ) GO  TO  410 
C 

C  CALCULATE  THE  DISTANCE  TO  THE  NEAREST  DXTRAN  SPHERE, 

DXL. 
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DXL=HUGE 

DO  80  1=1 ,NDX ( IPT) 

IF(IDX.EQ.I)GO  TO  80 
F=DXX ( IPT, 1 t I >-XXX 
G=DXX( IPT, 2, I ) -YYY 
H=DXX(IPT,3,I)-ZZZ 
Q=F*UUU+G*VW+H*WWW 
C=MIN(MAX  <  ZERO,Q) ,DXL) 

IF (  (F-UUU*C)**2+(G-VW*C)**2+<H-WWW#C>**2.LT.DXX(  IPT, 5,  I  )  > 

1 

DXL=MIN<  DXL , Q-SQRT (MAX <  ZERO , Q**2+DXX  < IPT , 5 , I ) -F**2-G**2-H**2 ) 
)  ) 

80  CONTINUE 
C 

C  CALCULATE  THE  DISTANCE  TO  TIME  CUTOFF,  DTC. 

DTC=VEL*  <  TCO (IPT) -TME  > 

C  CALCULATE  THE  CROSS  SECTIONS  IN  THIS  CELL. 

T0TM=0. 

PFP=0. 

STP=0. 

DEB=HUGE 

IF(MLL(LMLL+1 , ICL) .EQ.O>GO  TO  85 
IF ( IPT. EG. 1 ) CALL  ACETOT 
IF ( IPT. EQ. 2) CALL  PHOTQT 

SPECIAL  TREATMENT  FOR  MULT I GROUP  ELECTRONS. 

I F ( MC AL . EQ . O  >  GO  TO  83 
PFP=10.*PFP 
IF (STP.EQ.O. )G0  TO  85 
M=JXS< 1 ,MGEG8T ( 1 ) )+JGP-l 

I F ( MCAL . EO .  1  ) DEB=  ( ERG- YSS ( M )  ■*- . 5* YSS ( M+JGM (  1  )  )  )/STP 
IF ( MCAL .EQ .2 )DE8= ( YSS ( M ) + . 5* YSS ( M+JGM ( 1 ) >-ERG>/STP 

CALCULATE  THE  MEAN  FREE  PATH,  GS,  AND  ITS 
IPROCAL,  OPL. 

83  GS=0. 

PMF=HUGE 

QPL=  (  TOTM+PFP  >  #RHO  <  LRHO*  I  CL  > 

PLE=QPL 

IF ( PLE . EQ . 0 . ) GO  TO  160 
PFP=PFP/  ( TOTM+PFP ) 

IF (QAX (LOAX+IPT , I CL) .NE.O. )CALL  EXTRAN 
IFIOPL.LE.O. )G0  TO  160 
GS=1./QPL 

DECIDE  WHETHER  TO  FORCE  A  COLLISION. 

IF<FGR(LFOR+IPT, ICL) .EQ.0)G0  TO  90 
CALL  FORCOL 
IF(NTER.NE.O)GO  TO  310 
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GO  TO  160 


SAMPLE  THE  DISTANCE  TO  COLLISION,  PMF ,  NORMALLY, 
move  st  #  90  to  line  below 
90  CONTINUE 

qzridq=RANG ( > 
qzridql=-LOG<qzridq ) 

PMF=qzridql*GS 


THIS  PORTION  OF  CODE  CORRECTS  FOR  A  VARIABLE  DENSITY 
EXPONENTIAL  ATMOSPHERE 


CALCULATE  A  NEW  DISTANCE  TO  COLLISION,  PMF,  BASED  ON  A 
VARIABLE  DENSITY  ATMOSPHERE. 

NINP=0 

CALL  EQDIST (PMF,  D,  211,  WWW,  NINP) 

PARTICLE  EXITS  TOP  OR  BOTTOM  OF  ATMOSPHERE  IF  PMF=-1 
IF  (PMF.EQ.-l >PMF=HUGE 

RE-CALCULATE  THE  MEAN  FREE  PATH,  GS,  AND  ITS 
RECIPROCAL,  QPL. 

GS=PMF/qzridq 1 

QPL=1./GS 

PLE-QPL 


TALLY  THE  TRACK  LENGTH  IN  THE  CELL. 

160  D=MIN(PMF , DLS , DXL , DTC , DEB) 

I F ( NSR . NE .71) GO  TO  180 

I F  ( I PT .  NE .  1  .  OR .  LFCL  ( LLFC-*  I  CL  )  .  EO  •  0 )  GO  TO  180 
FM=0. 

DO  170  M=MLL<LMLLM  ,  I  CL)  ,MLL  (LMLL+2 , 1  CL) 

170  FM=FM+RTCR( 10,LME( 1 ,M> > #RTCR ( 8 , LME < 1 ,M>  >*FME(M> 
SUM><  <  3  >  =SUHK  <  3  >  +FM t  D t  WGT  *  RHO  ( LRHO+ 1  CL  ) 

180  L=LOCCT (LLCT+IPT , ICL) 

I F ( L . NE . 0 > CALL  TALLY(L,D> 

DO  185  1=0, LEV- l 

L=LOCCT (LLCT+IPT, INT (UDT (7,1 ) ) ) 

185  I F ( L . NE . 0 ) CALL  TALLY ( L , D ) 

JSU=JAP 


INCREMENT  THE  SUMMARY  ACCOUNTS. 

dt=d/vel 

PAC  ( LPAC-*  I  PT ,  5 ,  I  CL  )  =P  AC  ( LPAC-*  I  PT ,  5 ,  I  CL  >  *-WGT  *  DT  «ERG 
PAC  ( LPAC+ 1  PT  ,  6 , 1  CL  )  =PAC  C  LPAC*- 1  PT  ,  6  ,  I  CL  )  ♦WGTiDtERG 
PAC ( LPAC+ 1 PT , 7 , 1 CL  > -r AC ( LPAC+ 1 PT , 7 , I CL ) *D 
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I F ( PLE . NE . 0 . )PAC(LPAC+IPT ,8, ICL)*PAC(LPAC+IPT,B, ICL>+WGT*D/PL 
E 

PAC ( LPAC+ I PT , 9 , ICL  > -PAC ( LPAC+ 1 PT , 9 , 1 CL ) +WGT*DT 
PAC<LPAC+IPT,10, ICL)=PAC(LPAC+IPT, 10, ICL>+WGT*D 
C 

C  UPDATE  THE  PARTICLE  TO  THE  SURFACE,  COLLISION,  OR 

TERMINATION. 

C  BANKED  UNCOLLIDED  PART  COMES  BACK  HERE. 

190  TME=TME+DT 

XXX=XXX+UUU*D 

yyy-yyy+wv*d 

ZZZ=ZZZ+WWW*D 
DO  195  L=0, LEV-1 
UDT ( 1 , L ) =UDT ( 1 , L  >  +UDT ( 4 , L ) *D 
UDT ( 2 . L ) =UDT (2,L) +UDT ( 5 , L ) *  D 
195  UDT ( 3 ,L ) =UDT ( 3 , L ) +UDT ( 6 , L  >  *D 
C 

C  SPECIAL  TREATMENT  FOR  MULT I GROUP  ELECTRONS. 

I F ( STP . EQ . 0 . ) GO  TO  197 
T1=D*STP 

IF ( MCAL . EQ .2 ) T 1 =-T 1 
ERG=ERG-T 1 

PAX ( I PI ,6,2) =PAX ( IPT,6,2)+T1 *WGT 
IF(ERG.LE.ECF(1 ) .EQV.MCAL.LT. 2)G0  TO  270 
RM=YSS ( J  XS ( 1 , MGEGBT ( 1 ) >+JGP-l+2*JGM< 1 ) ) 

VEL=SLI TE4SQRT ( ERG* ( ERG+2 . *RM) ) / ( ERG+RM ) 

197  EG0=ERG 
C 

C  SCORE  FLUX  IN  CELL  FOR  MULT I GROUP  WEIGHT-WINDOW 

GENERATION. 

IF ( ICW . NE . 0 ) FLX ( LFLX ( IPT ) ♦MX A*  <  JGP-1 ) +  ICL  >  = 

1  FLX  <LFLX  (  IPT )  -t-MXA*  (JGP-1  >+ICL )  ♦D*WGT 
I F ( D . EQ . DTC ) GO  TO  2B0 
IFtD.EQ.DXLlGO  TO  300 
C 

C  PROCESS  DXTRAN  PARTICLE  AS  IT  LEAVES  ITS  SPHERE. 

IF ( I DX . EQ .0) GO  TO  200 

IF ( (XXX-DXX( IPT, 1 , IDX) ) *  *2-M  YYY-DX  X( IPT, 2, IDX) )**2+ 

1  <  ZZZ-DXX ( IPT, 3 ,IDX) )##2.LT.DXX(IPT,5,IDX) )G0  TO  200 
I  DX=0 

I F ( WWP < I PT ,4 ) . NE . O . ) GO  TO  200 
IF(WGT*FIM1  .GT.FIS*WCS2{  IPTMGO  TO  200 
T1=WCS1 ( IPT ) *FIS/F  ?M1 
IF ( T 1 . EQ.O. )G0  TO  200 
IF(WGT.LT.T1 *RANO{ ) )G0  TO  290 

PWB ( LPW8-*- IPT,  10,  ICL  >=PWB ( LPWB+IPT  ,10,  ICLJ+Tl-WGT 
PAX ( IPT ,2,6) =PAX ( IPT ,2,6) +T l-WGT 
PAXt IPT,3,6>=PAX<  IPT, 3 , 6  >■*•  <  Tl  -WGT  )  *ERG 
WGT=T1 
C 

C  ADJUST  THE  WEIGHT  FOR  EXPONENTIAL  TRANSFORMATION. 

200  IF 'QfiX (LQAX+IPT , ICL > .EQ .0. )G0  TO  210 
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T1=WGT 

WGT=WGT*EXP  < ( QPL-PLE ) *D ) 

I F ( PMF . LT . DLS >  WGT=WGT*PLE*GS 

PWB(LPWB+IPT,12, ICL)=PWB<LPWB+IPT, 12, ICL)-(Tl-WGT) 
1=2 

I F ( T 1 . GT . WGT ) I =5 

PAX  < I PT , I , 1 1 ) =PAX < I PT , I , 1 1 > +ABS ( T 1 -WGT ) 

PAX ( IPT , 1+1 , 11 )=PAX< IPT, 1+1 , 11)+ABS(T1-WGT)*ERG 
C 

C  PROCESS  THE  PARTICLE  THRU  THE  CELL  BOUNDARY  IF  NO 

COLLISION. 

210  I F  (  D . NE . DLS  >  GO  TO  220 
CALL  SURFAC 

I F ( KRFLG . NE . 0 ) CALL  EVENTP(3> 

IF (KDB . NE .0 > GO  TO  410 
I F ( NTER . NE . 0 ) GO  TO  310 
GO  TO  70 
C 

C  CALCULATE  EVERYTHING  ABOUT  THE  COLLISION. 

220  PAC ' LPAC+IPT , 3 , ICL>=PAC(LPAC+IPT,3, ICL5+1 . 

PAC ( LPAC+IPT , 4 ,  ICL)=PAC(LPAC+IPT,4,  ICD+WGT 

NCH ( IPT ) =NCH  < IPT  >  + 1 

NCP=NCP+1 

IF(NCH< IPT) ,EQ.DBCN(9) >G0  TO  260 
230  I F ( I PT . EQ . 1 > CALL  COLIDN 
I F ( I PT . EQ . 2 ) CALL  COLIDP 
I F ( KDB . NE . 0 ) GO  TO  410 
I F  <  NTER . NE . 0 ) GO  TO  310 


C  TALLY  DETECTORS  AND  CREATE  DXTRAN  PARTICLES. 

I F ( NDET (IPT) . NE . 0 ) CALL  TALLYD 
I F ( KDB . NE . 0  >  GO  TO  410 
IF (NDX ( IPT) .EQ.O)GO  TO  255 
IF <NDX< IPT) .GT. 1 .OR. IDX . EQ . 0 > CALL  DXTRAN 
I F ( KDB . NE . 0 ) GO  TO  410 
C 

C  PLAY  THE  WEIGHT-WINDOW  AND  ENERGY-SPLITTING  GAMES. 

255  IF(ABS(WWP( IPT, 4) > .EQ. 1 . .AND. I DX . EQ . 0 ) CALL  WTWNDO(WW) 
IF (NTER .NE . 0  >G0  TO  310 
I F ( ESPL ( I PT i 1 ) . NE . 0 . ) CALL  ERGIMP 
I F ( NTER . NE . 0 ) GO  TO  310 
GO  TO  60 
C 

C  DEBUG  FEATURE:  COLLISION  LOOP  BREAKPOINT. 

260  GO  TO  230 
C 

***************  PROCESS  TERMINATED  PARTICLES. 
**************** 

* 

270  NTER=4 
GO  TO  310 
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280  NTER=13 
GO  TO  310 
290  NTER=6 

GO  TO  310 
300  NTER=10 
C 

C  INCREMENT  PARTICLE  STATISTICS  FOR  TERMINATION  TYPE 

NTER. 

310  I F ( KRFLG . NE . 0 ) CALL  EOENTP  <  5 ) 

I F  <  I GWW . NE .  0  .  AND .  <  NTER . LT  . 6 . GR . NTER .  GT  . 9 )  )CALL 
WGTWWG ( 1 , WGT ) 

j=jrwb(nter; 

I F ( J . NE . 0 ) PWB ( LPWB+ I PT , J , I CL ) =PWB ( LPWB+IPT , J , I CL  > -WGT 
IF  <NTER.EQ. 1 >TMAV( IPT, 1 )=TMAV( I PT , 1 ) +TM£*WGT 
IF\NTER.EQ.3)TMAV( IPT,2)=TMAV( IPT , 2 ) +TME*WGT 
TMAV < I PT , 3  >  =TMAV < I PT , 3  >  +TME* WGT 
PAX ( IPT ,  4 , NTER ) =PAX  < IPT , 4 , NTER )  +  1 . 

PAX ( I PT  >  5 , NTER ) =PAX ( I PT , 5 , NTER ) +WGT 
PAX  < I PT , 6 , NTER  >  =PAX ( I PT , 6 , NTER ) +WGT*£RG 
IF<NSR.NE.71 >60  TO  320 
IF ( IPT . NE . 1 > GO  TO  320 

I F  <  NTER . NE . 1 . AND . NTER . NE . 4 . AND . NTER . NE . 1 3 ) GO  TO  320 
RLT ( 1 ) =RLT ( 1 >+WGT*TME 
RLT ( 2 ) =RLT ( 2  >  +WGT*TME 
320  NT£R=0 
C 

C  GET  THE  NEXT  PARTICLE  FROM  THE  BANK,  IF  THERE  ARE 

ANY. 

I F ( NBNK . EQ . 0  >  GO  TO  390 
DO  330  1=1 ,MXA 
330  I PAC2 ( LPC2+ I ) =0 
CALL  BANK I T  C  2 ) 

100  IF( KFL . EG. 0 > GO  TO  370 
IF (KFL.EQ. 1 )GG  TO  350 
DO  340  1=1, MXA 

340  IF ( IFL(LIFL+I ) . GE . NODE  > I FL < L I FL+ 1 >=0 
IF ( KFL.EQ , 2 >G0  TO  370 
350  DO  360  J=i,MXJ 

360  I F (JFL ( LJFL+J  > . GE . NODE > JFL < L JFL+J  > =0 
370  IF ( KRFLG. NE.O) CALL  EVENTP(2> 

PAC(LPAC+IPT,  2,  I  CL  )  =PAC  (LPAC-*-IPT  ,2, 1CU+1  . 

I PAC2  <  LPC2* I  CL  >  =  1 
IF{NPA.LT.0)G0  TO  380 
C 

C  PROCESS  PARTICLE  FROM  THE  SURFACE  SOURCE. 

I F  <  J  SU . GE . 0 ) GO  TO  60 

IF(ND£T (IPT) .EQ.O. AND.NDX < IPT  > . EQ . 0 ) GO  TO  375 
1PSC=12 
SWTM=WGT 
CALL  STARTP 
375  JSU=ABS(  JS'J) 

IF(WC1 (IPT) .6T.0. >G0  TO  60 
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WCS1 ( IPT)=-WC1 < IPT>*WGT 
WCS2  < I PT ) =-WC2  < I PT  >  * WST 
GO  TO  60 
C 

C  SHORT  LOOP  IF  PARTICLE  IS  UNCOLLIDED  PART  OF  FORCED 

C  COLLISION. 

380  D=MIN<DLS,DXL,DTC> 

DT=D/VEL 
PMF=HUGE 
JSU=NPA+1 000000 
GO  TO  190 
C 

C  THE  HISTORY  IS  COMPLETE. 

C  ADD  THE  TALLY  DATA  OF  THIS  HISTORY  TO  THE  TOTAL 

TALLY  DATA. 

390  NCT ( 1 ) =NCT ( 1 ) +NCH ( 1 ) 

NCT ( 2 ) =NCT ( 2 ) +NCH ( 2 ) 

I F ( MODE . EQ . 2  >  GO  TO  400 
SMUL ( 2 ) =SMUL  <  2 ) +SMUL ( 1 ) 

SMUL  <  3 ) =SMUL ( 3  >  +SMUL ( 1 ) *  *  2 

400  I F ( NT AL . GT . 0  >  CALL  TALSHF 

IF (NTAL .GT.O . AND .MODINPS, NPD ) .EQ.O )CALL  ADDTFC 
RANK=RANI 
RANL=RANJ 
RETURN 
C 

*******************  PROCESS  LOST  PARTICLE. 
******************* 

* 


410  I F ( KDB . GE . 1 1 >  GO  TO  450 
C 

C  CLEAR  THE  FIRST  TALLY  BLOCK  AND  READ  VARIABLE  COMMON 

BACKUP. 

DO  420  1=1, MXF 
420  TAL(LTAL+I )=0. 

I F ( NRSW . NE . 0 ) CALL  SUFWRT ( 0 , ZERO) 

DO  425  1=1 , NVARCM 
425  6VARCM ( I ) =GVBU ( I ) 

DO  430  1=1 , LVARCM 
430  JVARCMC I )=JVBU< I ) 

C 

C  RERUN  HISTORY  WITH  FULL  SURFACE  SENSE  CHECK  AND 

EVENT 

C  PRINTING. 

IF (KRFLG.EQ . 2 .OR .NERR.GE .LOST ( 2 ) ) GO  TO  440 
KRFLG=2 
NT  1 1=0 
GO  TO  40 
C 

C  RETURN  TO  PRINT  DEBUG  INFORMATION  AND  START  A  NEW 

HISTORY. 

440  NERR=NERR*1 
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on  n  n  n  n  n 


T  F  <  NERR . LE . LOST  <  2 ) >K0V=1 
RETURN 
C 

******************  TERMINATE  LONG  HISTORY. 
******************* 

* 

450  NST =NST +256 

I F ( NRSW . NE . 0 ) CALL  SUFWRT <  0 , ZERO  > 

DO  460  1=1 , NVARCM 
460  GVARCM ( I ) =GVBU ( I ) 

DO  470  1=1 , LVARCM 
470  JVARCMd  >=JVBU<  I  > 

RANI=RANK 
RANJ=RANL 
IF ( NSR .NE .6)  then 
RETURN 
end  i  -f 

BACKSPACE  I USR 
4B0  BACKSPACE  I USR 
BACKSPACE  I USR 
NRRS=NRRS-i 
READ( I USR) A 

I F ( NPSR . EQ . ABS ( A  >  >  GO  TO  4B0 

NQSS=NQSS-1 

RETURN 

END 


SUBROUTINE  TRANSM ( DD , ST > 

C  CALCULATE  THE  ATTENUATION  AMFP  OYER  THE  DISTANCE  DD 

FROM 

XXX  ,  YYY  ,  ZZZ  IN  THE  DIRECTION  UUU,VVUfWWW. 

ST  IS  THE  RUSSIAN  ROULETTE  LEVEL, 
include  ’CM.inc' 
include  'RC.inc' 

SD=0. 

AMFP=0. 

IF<LCA<LLCA+ICL> .LT.OJCALL  CHKCEL (I CL , 3 , J ) 

FT=ST 

IF ( FT . NE . 0. >TT=-LOG<FT> 

DO  RUSSIAN  ROULETTE  ON  SMALL  SCORES. 

10  IFfFT.EQ.G. >G0  TO  20 
IF(AMFP.LT.TT)GO  TO  20 
T=EXPC-AMFP> 

IF<FT*RAN6( ) .GT.T)GO  TO  30 
WGT=WGT*FT/T 
FT=T 

tt=amfp 

CALCULATE  THE  ATTENUATION  FOR  THIS  SECTION  OF  THE 


(J*****-**+*UU  L)  CJ  (J  -»•»*•»*  (J  «  m  *  u  U 


TRACK . 

20  CALL  TRACK < I CL ) 

IF(KDB.NE.O)  then 
RETURN 
end  i-f 
T0TM=0. 

I F ( MLL ( LKLL+ 1 , 1 CL  > . EQ , 0  >  GO  TO  25 
I F  < I PT . E  J  .  1 >  CALL  ACETOT 
IF ( IPT.EQ.2)CALL  PHOTOT 
25  PLE=TOTM*RHO<LRHO+ICL> 

C 

D=M I N ( DLS , DD-SD ) 


THIS  PORTION  OF  CODE  CORRECTS  FOR  A  VARIABLE  DENSITY 
ATMOSPHERE 


CALCULATE  A  NEW  MACROSCOPIC  CROSS  SECTION,  PLE,  BASED  ON 
A  VARIABLE  DENSITY  ATMOSPHERE. 

NINP=1 

TAKE  RECIPROCAL  OF  MACROSCOPIC  CROSS  SECTION 
XMFP=1 ./PLE 

CALL  EQDIST ( XMFP,  D,  ZZZ,  WWW,  NINP) 

NINP=0 

PLE=1./XMFP 


CALCULATE  A  NEW  NUMBER  OF  MEAN  FREE  PATHS,  AMFP . 
AMFP=AMFP+PLE*D 


I F  <  AMFP , GT . 80 . ) GO  TO  30 
SD=SD+DLS 

UPDATE  THE  LOCATION  AND  THE  NEW  CELL. 
XXXsXXX+UUUKD 
YYY=YYY+VVV#D 
ZZZ=ZZZ+WWW*D 
DO  27  L=0,L£V-1 
UDT  {  1  ,  L  )  =UDT  <  1 ,  L )  -HJDT  <  4  ,  L )  *  D 
UDT <  2 , L ) =UDT <  2 , L  >  +UDT  C  5 , L  >  *  D 
27  UDT ( 3 ,L >  =UDT ( 3,L >+UDT ( 6, L) $D 
TME=TME+D/VEL 
IF(SD.GE.DD)  then 
RETURN 
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n  n 


* 


end  if 
JSU=JAP 

CALL  NEWCEL(ZERO) 

IF(KDB.NE.O)  then 
RETURN 
end  if 
ICL=I AP 

IF ( F I M  <  LF I M+ 1 PT , I AP ) . NE . 0 . ) GO  TO  10 
CALL  BEY0ND(5> 

RETURN  WITH  ZERO  WEIGHT  FOR  ZERO  IMP.  OR  IF  SCORE  IS 

REJECTED. 

30  WGT=0. 

RETURN 

END 
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